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EXECUTIVE  SUMMARY 


This  report  presents  an  approacn  to  describing  tne  performance 
of  mobile,  maintenance  contact  teams  (C-teams).  Tnese  teams  support 
an  operational  system  within  a  specific  service  area.  In  the  interest 
of  generality,  the  supported  system  is  left  undefined.  However,  values 
selected  for  a  parametric  analysis  were  based  on  a  tactical  air  def¬ 
ense  system.  Tne  ultimate  objective  of  C-team  service  is  to  Keep 
tne  customer  system  in  a  high  state  of  operational  readiness.  Tnus, 
C-team  performance  is  measured  by  variables  such  as  average  number 
of  nonoperational  customer  units  and  average  time  a  failed  unit  stays 
nonopera tional .  Although  the  mathematical  model  of  the  system  in 
steady-state  is  fully  analytic,  it  is  not  an  expected-value  model. 

Tne  probability  distribution  for  number  of  nonopera tional  customer 
units  is  an  output  of  tne  implementing  computer  program.  Statistics 
summarizing  the  performance  of-  and  utilization  of  C-teams  are  also 
outputs . 

A  parametric  analysis  was  performed  to  examine  the  sensitivity 
of  C-team  performance  to  parameters  such  as:  (a)  intercustomer 
speed  and  variation  in  speed,  (b)  density  of  C-teams,  (c)  dedicated 
versus  nondedicated  C-teams,  (d)  mean  time  per  unit  between  service 
requests,  (ej  mean  diagnostic  and  repair  time,  and  (f)  form  of  the 
probability  distribution  for  total  service  time.  Some  conclusions 
are  drawn  from  tnese  results,  which  may  be  applicable  to  an  oper¬ 
ational  system  of  interest.  For  example,  when  each  of  3  C-teams  is 
dedicated  to  1/d  of  the  customer  population,  a  substantial  perform¬ 
ance  penalty  is  paid  relative  to  an  arrangement  of  nondedicated 
C-teams,  which  serve  the  entire  population. 

For  the  interested  analyst,  the  report  sketches  the  method  of 
derivation  of  formulas.  Computer  source  code  for  the  implementing 
programs  is  included  with  extensive  comments. 
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MEMORANDUM  REPORT 

SUBJECT:  Analytic  Simulation  of  the  Performance  of  Mobile 
Maintenance  Contact  Teams 


1 . References 


a.  Memorandum  Report  No  DRSMC/SA/MR-4 ,  (AD  A136580),  HQ,  AMCCOM, 
Nov  8j,  title:  Models  of  a  Service  System  for  Production  Machine 
Maintenance. 

b.  Technical  Report,  Society  for  Industrial  and  Applied  Math, 
1979,  title:  LINPAC  Users'  Guide. 

c.  Textbook  by  Donald  Gross  and  Carl  Harris,  J.  Wiley,  c.  1974, 
title:  Fundamentals  of  Queueing  Theory. 

2.  Background 


A  maintenance  contact  team  (C-team)  is  an  adjunct  to  Organizational- 
and/or  DS-  maintenance.  Each  C-team  is  considered  a  mobile  service 
system  which  operates  within  a  prescribed  service  area.  The  C-team 
concept  treated  here  is  a  derivative  of  the  radar/elec tronics  van 
supporting  the  Vulcan  Air  Defense  System,  wnich  was  developed  in 
the  1960's.  A  C-team  is  dispatched  so  that  it  can  move  to  the 
next  customer  immediately  upon  completing  service  on  the  current 
customer.  In  addition  to  the  time  to  diagnose  and  repair  a  fault, 
the  total  time  a  C-team  "serves"  a  customer  includes  the  time  to 
travel  to  the  customer,  following  a  commitment  to  tnat  customer. 
Intercustoraer  distance  and  speed  are  variables  which  affect  tne 
service  performance.  Variations  in  these  random  variables  as 
well  as  the  conventional  variation  in  maintenance  service  time 
are  factors  considered  in  modeling  this  service  system.  Per¬ 
formance  of  this  service  system  is  measured  by  the  operational 
availability  of  the  system  being  supported.  To  be  general, 
the  latter  system  is  not  specified  here. 

3.  Objectives 


An  objective  of  this  MFR  is  to  describe  an  analytic  --  as 
opposed  to  Monte-Carlo  —  model  wnich  simulates  the  stochastic 
steady  state  of  a  service  system  consisting  of  (possibly)  multiple 
C-teams  serving  a  finite  population  of  customers  within  a  service 
area.  A  secondary  objective  is  to  present  some  parametric  results 
wnich  may  be  helpful  in  describing  tne  adequacy  of  C-team  performance. 
Parameter  values  used  in  this  analysis  are  based  on  a  range  of 
values  considered  applicable  to  the  SGT  York  air  defense  system. 

Tne  relative  importance  of  certain  parameters  and  insensitivity 
of  others  are  displayed. 


4.  Model  Assumptions 


The  customers  are  located  witn  uniform  probability  density  in 
a  rectangular  service  area.  Contact  teams  are  dispatched  from 
one  customer  to  the  next  without  leaving  the  service  area, 
m/aen  a  o-team  completes  service  on  one  customer,  it  either  (a) 
immediately  travels  to  the  next  customer,  if  a  customer  is  waiting, 
or  ( b)  waits  at  current  position  until  dispatched,  if  the  queue 
of  customers  is  empty.  Tne  path  traveled  by  C-teams  between 
customers  follows  a  series  of  segments,  each  of  which  is  parallel 
to  one  of  the  sides  of  the  service  area.  (The  last  assumption  rep¬ 
resents  use  of  road  networks  and  the  need  to  avoid  obstructions. ) 
The  intercustomer  speed  varies  randomly  over  occasions  witn  a 
uniform  probability  density  between  prescribed  limits.  Tne 
next  customer  to  be  served  is  selected  by  means  of  one  of  two 
service  disciplines  --  FIFO  or  SDST.  Using  FIFO,  the  next 
customer  is  the  one  whose  request  came  first.  Using  SDST,  the 
next  customer  is  the  one  with  the  shortest  distance  from  the  team's 
current  position.  All  prametric  results  shown  here  use  FIFO 
service  discipline.  (When  the  average  queue  length  is  small,  there 
is  no  practical  difference  in  results  from  the  two  disciplines.) 

The  requests  for  C-team  service  are  assumed  to  occur  at  random  in 
a  conditionally  Poisson  manner.  The  rate  parameter  per  operational 
unit  is  a  constant.  Tnus,  arrival  rate  of  requests  is  proportional 
to  the  number  of  operational  units,  i.e.,  units  not  down  for  C-team 
main tenance .  Disabling  faults  other  than  those  repairable  by 
C-teams  are  not  treated  explicitly  in  this  model.  The  diagnostic 
and  repair  function  is  always  assumed  to  restore  a  failed  unit  to 
an  operational  status.  The  average  time  to  do  this  is  the  MTTR. 

The  nominal  values  of  the  variables  used  in  this  study  are  given 
in  Table  1.  Parametric  excursions  were  made  for  selected  variables. 

TABLE  1 

VARIABLES  USED  IN  THE  CONTACT-TEAM  PARAMETRIC  ANALYSIS 
Standard  Parameter  Values 


Frontal  Width  of  Service  Area  _  20  km  * 

Cross -front  Dimension  of  Area  4  km 

Intercustomer  Average  Speed  20  km/nr 

Range  of  Uniformly  Dist'd  Speed”-  5  km/hr 
Avg  Time  Between  Service  Requests  100  hr/unit 
Avg  Diagnotic  and  Repair  Time  hr 

Population  Dens  of  Fire  Units  to  36 

Number  of  Supporting  C-Teams  1  to  3 


Front  of  service  area  represents  a  Division  Front. 


b.  Methodology 


Every  aspect  of  this  problem  has  been  approached  analytically, 
and  has  been  verified  by  Monte-Carlo  methods.  This  includes  the 
model  of  C-team  movement  between  customers  as  well  as  the  stocnastic 
service  or  queueing  model.  One  advantage  of  a  purely  analytic  ap¬ 
proach  is  that  parametric  analyses  can  be  performed  with  assurance 
that  tne  observed  differences  in  output  statistics  are  in  no  way 
affected  by  Monte-Carlo  sampling  variation.  Tnis  aspect  is  parti¬ 
cularly  important  when  the  natural  variation  in  subject  variables 
is  relatively  large  or  wnen  differential  parametric  effects  are 
small.  Altno  some  of  the  model  assumptions  may  seem  somewhat  restr¬ 
ictive  —  such  as  a  uniform  distribution  of  customers  within  the 
service  area,  these  have,  in  fact,  been  found  not  to  be  so.  Some 
alternatives,  which  are  not  analytically  tractable,  have  been  ex¬ 
amined  via  Monte  Carlo.  These  auxiliary  studies  are  not  reported 
here. 

6.  Tne  service  performed  by  a  C-team  has  two  stages  --  travel  to 
a  customer  and  repair  of  a  fault.  Therefore,  it  is  appropriate  to 
consider  as  a  model  of  the  service  system  a  Markov  model  having 
two  stages  of  service  and  serving  a  finite  customer  population. 

The  total  time  to  "serve"  a  customer  is  considered  a  gamma(2) 
random  variable.  This  model  is  equivalent  to  two  identical  stages 
of  exponential  service  in  tandem.  Alternative  models  of  service 
were  also  used.  Tnese  are  discussed  under  "Results".  A  multi- 
server  analytic  model  of  this  sort  was  used  earlier  (Ref  a.)  to 
represent  the  service  system  for  production  machine  maintenance. 
Numerical  methods  are  used  to  solve  the  steady-state  equations, 
which  derive  from  the  Kolmorogov  diff erential-difference  equations 
for  tne  state  probability  vector.  In  the  present  application,  it 
was  considered  necessary  to  change  the  implementing  computer 
program  in  order  to  reduce  execution  time  for  cases  naving  a  large 
customer  population.  Program  changes  responsible  for  a  great  red¬ 
uction  in  CPU  time  are:  (a)  use  of  a  better  method  for  solving  the 
state  equations  and  (b)  reduction  in  the  number  of  Markov  states. 

In  the  computer  program,  found  in  Annex  B,  the  reduced  state- 
transition  matrix,  derived  from  the  Kolmogorov  equations,  is,  first, 
factored  (into  upper  and  lower  triangular  matrices),  and  then  invert¬ 
ed  in  place.  Two  LINPAC  routines  (Ref  D. )  are  used  for  this  purpose. 
This  approach  is  about  2. 5  times  faster  than  the  gaussian-reduction 
method  of  Ref  a.  for  large  (  g.t.  TOO  X  100)  matrices.  The  dimension 
of  the  state  probability  vector  is  reduced  by  eliminating  states 
with  negligible  probability  of  being  occupied.  The  choice  of  states 
to  eliminate  is  based  on  the  closed-form  solution  to  a  similar 
queueing  problem.  The  similar  problem  involves  a  finite  customer 
population  with  conditionally  Poisson  arrivals  and  exponential 
services,  whose  steady-state  solution  is  found  in  Reference  c. 

The  deleted  states  are  those  associated  with  a  number  in  the  service 
system  greater  than  N,  where  the  probability  that  the  system  number 
exceeds  N  is  smaller  than  some  small  value  (acceptable  error). 

An  error  probability  of  1/10,000  generally  reduces  the  number  of 
states  by  more  than  a  factor  of  1/2  for  the  cases  treated  here. 

This  reduction  is  quite  significant,  since  CPU  time  is  nearly 
proportional  to  the  number  of  states  squared. 
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7.  Survey  of  Results 


Several  formulas  of  interest  were  derived  for  this  model. 

Tne  system  performance  requires  inputs  to  the  queueing  model 
sucn  as:  (a)  average  time  a  C-team  spends  in  travel  between 
customers  and  (b)  average  maintenance  service  time  per  customer. 

Tbe  first  variable  depends  upon  the  mean  and  variance  of  the 
intercustomer  distance  and  the  mean  and  variance  of  the  inter¬ 
customer  speed.  Formulas  were  derived  for  the  probability 
distribution  of  intercustomer  distance  and  for  the  mean  and 
variance  of  this  distribution.  Since  the  intercustomer  speed 
is  considered  a  uniform  random  variable  over  the  range  (si,  s2), 
tbe  mean  and  variance  of  the  speed  are  tne  standard  results: 

(si  +  s2)/2  and  (s2  -  s11**2/12,  respectively.  An  approximation 
for  the  mean  and  variance  of  the  ratio  of  two  random  variables 
can  be  used  to  obtain  the  mean  and  variance  of  the  intercustomer 
travel  time  —  =  distance/speed  —  in  terms  of  the  means  and 
variances  of  distance  and  speed.  The  accuracy  of  the  approx¬ 
imation  depends  somewhat  upon  the  skewness  of  the  distributions 
of  distance  and  speed.  The  approximation  was  checked  against 
an  exact  metnod  for  a  practical  range  of  values  of  the  distance 
and  speed  parameters  and  found  to  be  quite  adequate.  The  approx¬ 
imation  can  be  given,  simply,  as  follows.  Let  E(d)  and  V(d) 
be  the  mean  and  variance  of  the  intercustomer  distance.  Similarly, 
let  E(s)  and  V(s)  be  the  mean  and  variance  of  the  speed.  Tnen,  the 
mean  and  variance  of  the  travel  time  are  given,  respectively,  by 

3 

=  E(d)/E(s)  +  E(d) V (s) /( E(s1 ) 

2  2  4  2  i 

=  V(d)/(E(s) )  +  V(s) (E(d) )  /(E(s) )  +  (E(d)V(s))  /(E(s)) 

4 

+  V(d)V(s)/(E(s) )  . 

Equations  for  the  mean  and  variance  of  the  intercustomer  distance 
for  a  FIFO  service  discipline,  are  given  in  terms  of  the  dimensions 
of  the  rectangular  service  area,  say,  a  by  b: 

E(d)  =  (a  +  b)/j 

and 

2  2 

V(d)  =  (a  +  b  J/18  . 

Equations  for  the  cumulative  distribution  function  (C.D.F..)  of 
intercustomer  distance  are  presented  in  Annex  A.  The  method  of 
derivation  is  outlined  in  this  annex.  Annex  B  contains  a  computer 
source  listing  of  the  program  which  implements  all  the  auxiliary 
equations,  such  as  tnose  above,  as  well  as  the  equations  of  the 
Markov  queueing  model.  An  example  of  the  C.D.F.  of  intercustomer 
distance  is  given  in  Figure  1  using  the  nominal  dimensions  of 
the  service  area.  The  C.D.F.  of  intercu3tomer  distance  using  a 
SDST  service  discipline,  given  n  customers  are  queued  (n  g.t.  0), 
is  also  shown  in  Figure  1.  Notationally ,  let  the  latter  C.D.F. 


EU) 

and 

V  ( t) 
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oe  G(n, d)  and  the  corresponding  C.D.F.  using  FIFO  be  F(d),  for 
intercustomer  distance  d.  Then,  these  probaoility  distributions 
are  related  as  follows: 


n 

G(n,d)  =  1  -  (1  -  F( d) ) 

8.  Tne  auxiliary  equations  yield  the  mean  total  time  in  which  a 
C-team  is  involved  with  a  customer.  Tnis  parameter  and  the  mean 
time  between  failures  of  operational  units  (MTBF)  are  inputs 

to  the  Markov  queueing  model  (Annex  B).  Tne  outputs  of  this 
model  are  typical  of  stochastic  service  models:  (a)  Average  and 
standard  deviation  of  numoer  of  individuals  in  the  service  system, 
i.e.,  either  waiting  (in  a  service  queue)  for  a  C-team  to  be  assigned 
or  being  "served”.  In  this  context  "being  served"  means  that  a  C-tea a 
is  either  in  transit  or  conducting  repairs,  whereas  the  expression 
"in  the  service  system",  used  in  queueing  theory,  means  that  an 
individual  operational  unit  has  experienced  a  failure  (of  the  type 
which  is  repairable  by  a  C-team),  wnich  has  not  been  repaired. 

In  this  context,  "in  tne  service  system"  just  means  "nonopera tional" . 

(b)  Avg  and  S.D.  of  number  of  individuals  in  the  service  queue.  In 
this  situation  there  is  only  one  service  queue,  since  failed  units 
request  service  from  a  single  source,  the  dispatcher. 

(c)  Avg  and  S.D.  of  number  of  busy  service  channels  (  or  busy  C-teams ) . 

(d)  Avg  and  S.D.  of  waiting  time  in  the  service  system,  i.e.,  the 
time  a  customer  is  "down"  with  a  failure,  repairable  by  a  C-team. 

(e)  Avg  and  S.D.  of  waiting  time  in  the  service  queue.  This  time 
interval  starts  with  a  request  for  service  and  ends  when  a  C-team  is 
committed  to  serve  the  requesting  unit. 

9.  Model  runs  were  made  under  two  contingencies  (or  options; 

for  a  battalion  population  of  jo  operational  units.  In  the  first 
option  the  3o  units  are  divided  into  batteries  of  12  units,  each 
of  wnich  has  a  single  C-team  dedicated  to  serving  that  battery. 

In  the  second  option  any  of  the  3  C-teams  can  serve  any  customer. 

Tne  latter  option  does  not  "fence  off"  requests  by  the  artifice 
of  restricting  tne  population  which  a  C-team  can  serve. 

10.  Using  a  (battery)  population  of  12  units,  the  effect  of 
in tercustomer  speed  on  maintenance  performance  is  shown  in 
Table  2.  The  average  speed  is  increased  from  10  to  40  km/hr, 
in  the  first  three  runs  —  columns  1  thru  3.  In  the  fourth 

run  the  range  of  speed  is  increased  while  preserving  the  average 
at  the  nominal  value,  20  km /hr.  Tnis  may  be  compared  with  the 
results  in  column  2  to  determine  the  effect  of  variability  in 
the  speed.  In  the  fifth  run,  the  density  of  C-teams  is  tripled 
to  shown  density  effects. 
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TABLE  2 


PARAMETRIC  EFFECTS  OF  CONTACT-TEAM  MOBILITY  AND  DENSITY  ON  SERVICE 


Attribute  Service  System  Parameters:  Teams  per  Btry,  Range  of  Speed 


ui  out:  - 

Performance 

i 

i  i 

i  + 

1  o 

1,  20+-5 

1,  40+-5 

1,  20+-1 0 

3,  20+ -5  km/n; 

Avg  Num  in  Sys* 

0.648 

0.546 

0T503 

0.552 

0.399 

Avg  Units  in  Q 

0.209 

0.15b 

0.153 

0.159 

0.001 

Avg  Wait  in  Sys 

5.712 

4.770 

4.375 

4.819 

3.439  Hrs 

Avg  Wait  in  Q 

1.845 

1.362 

1 .174 

1.385 

0.00o  Hrs 

Pr(Q  Wait[8  hr) 

0.933 

0.961 

0.970 

0.959 

0.999 

*  "Average  Customers  in  the  Service  System"  is  equivalent  to  the 
number  of  customer  units  "down"  for  service  by  a  contact  team. 

11.  Using  the  second  option,  a  series  of  runs  were  made  with 
MTbF  and  MTTR  as  parameters.  Tne  excusrions  in  tnese  parameters 
were  chosen  to  develop  the  locus  of  points,  in  the  space  of  these 
parameters,  on  which  the  average  number  of  nonopera tional  customer 
units  has  the  constant  value  —  0.5,  1.0,  1.5,  or  2.0.  This  family 
of  curves  is  displayed  in  Figure  2.  These  curves  are  isoavail¬ 
ability  curves  having  availabilities  of  98. 6%,  97.2%,  95.8%,  and 
94.4%  for  the  supported  system,  respectively.  An  alternative  way 
of  displaying  these  data  is  shown  in  Figure  3.  In  this  figure 
the  average  number  of  customer  units  "down"  is  snown  as  a  function 
of  MTBF,  with  MTTR  as  a  parameter.  Note  that  the  curve  for 
each  particular  MTTR  shows  a  "Knee"  at  whicn  the  rate  of  change 
of  average  nonopera tional  units  changes  remarKably. 

TABLE  3 

SERVICE  PERFORMANCE  FOR  A  SYSTEM  OF  3  C-TEAMS  SERVING  3b  CUSTOMERS 


!  Avg  Customers  in  Service  Sys  and  Avg  Wait  (Hrs)  in  Service  Sys 
MTBF  I  for  an  MTTR  (Hrs) : 

( Hrs)  !  2  i  3  |  4  j  5 


40“ 

!  2. 558 

3.059” 

60 

!  1  .<1*6 

2.6 02  ! 

75 

j 

1 

80 

r  1.088 

2*4^4  ! 

100 

£  0.862 

2.454  ! 

120 

t  0*7  16 

2.436  3 

125 

r 

! 

150 

f 

] 

160 

!  Uo38 

2.421  J 

200 

!  0*4  ju 

2.415 

250 

; 

300 

! 

350 

1 

400 

! 

! 

! 


2.342 

4.175  ! 

1 

1.736 

3.800  ! 

1 

!  2.468 
| 

1  .244 

3.580  ! 

| 

!  1.672 

1 

o.ybl 

3.500  ! 

i  1.292 

0.812 

3.464  ! 

| 

!  1.061 

| 

0.608 

3.433  ! 

!  0.788 

0.486 

3.4c:2  ! 

!  0.629 

0.405 

3.416  ! 
1 

0.524 

0.449 

! 

I 


! 

I 

5.512  ! 


] 

4.872  ! 

1 

1 

!  2.183 

1 

6.454 

4.653  ' 

1 

4.555  ! 

1 

!  1.324 

1 

5.729 

4.474  ! 

!  0.972 

5.561 

4.444  ! 

!  0.773 

6 . 486 

4.429  ! 

0 . 64  j 

5.464 

4.422  ! 

0.551 

6.4  jb 

1 

0.482 

5.429 
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Figure  2.  Customer  Isoavailability  Curves  ir>  the  Space  of  MTBF  and  MTTR 
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Figure  3.  Average  Nonoperational  Customer  Units  as  a  Function  of  MTBF  with  MTTR  as  a  Parameter 


12.  The  results  presented  above  were  obtained  for  one  statistical 
model  of  the  time  for  a  C-team  to  "serve"  a  customer,  i.e.,  to  reach 
a  customer  and  to  diagnose  and  repair  the  fault.  This  model 
uses  a  gamma  probability  distribution,  with  shape  parameter  of  2. 
This  is  denoted  as  the  gamma(2)  model.  Among  the  alternative 
models  one  might  reasonably  take  are:  (a)  the  service  activity 
time  is  exponentially  distributed,  or  (b)  the  travel  time  is 
exponentially  distributed  (as  a  1st  stage  of  service),  and  the 
actual  repair  time  is  exponentially  distributed,  as  a  2nd  stage 
of  service.  To  examine  the  sensitivity  of  the  queueing  results 
to  the  distributional  assumption,  comparison  runs  were  made  for 
the  6amma(2)  and  exponential  models.  A  gamma(2)  model  is  actually 
equivalent  to  two  tandem  exponential  models  having  identical  rate 
parameters.  Tnerefore,  the  gamma(2)  and  exponential  models 
produce  results  which  bound  results  of  model  (b),  when  all 
models  are  constrained  to  have  the  same  mean  activity  time.  The 
results  of  the  paired  comparisons  are  shown  in  Table  4.  Note 
that  when  the  C-teams  are  relatively  inactive,  the  model  alterna¬ 
tives  produce  results  whicn  are  quite  close,  particularly  with 
respect  to  the  average  number  of  units  in  the  service  system. 

Tne  values  of  parameters  which  are  treated  as  constants  in  Table  4 
are : 

Dimensions  of  the  Service  Area  20  km  X  4  km 

Intercustomer  Speed  20  +-  5  km/nr 

Customer  Population  Density  ^6 

Numoer  of  Supporting  C-Teams  3 

Tne  two  parameters  shown  in  Table  4  are  MTBF,  the  mean  time  (hours) 
between  service  requests  per  unit,  and  MTTR,  the  mean  diagnostic 
and  repair  time  (hours). 
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TABLE  4 


COMPARISON  OF  QUEUEING  STATISTICS  FOR  ALTERNATIVE 
DISTRIBUTIONS  OF  C-TEAM  TRAVEL-AND-SERVICE  TIME 


C.D.F.  Statistical  Attribute  of  C-Team  Performance 


Total 

Time 

Avg  No 
in  Sys 

S.D.  No 
in  Sys 

Avg  No 
in  Queue 

Avg  Busy 
Servers 

S.D.Busy 

Servers 

Avg-Sys 
Wait  Tm 

Avg-Queue 
Wait  Time 

MTBF  =  60 
MTTR  =  1 

Gamma (2 ) 
Exponen 

0.8400 

0.8433 

0.9213 

0.9347 

0.0147 

0.0181 

0.8253 

0.825B 

0.8737 

0.8748 

1.433 

1.439 

0.025 

0.031 

MTBF  =  7b 
MTTR  =  1 

Gamma (2 ) 
Exponen 

0.6697 

0.6711 

0.8183 

0.8247 

0.0063 

0.0078 

0.6634 

0.8634 

0.7946 

0.7953 

1.422 

1.425 

0.013 

0.01b 

MTBF  =100 
MTTR  =  1 

Gamma (2 ) 
Exponen 

0.5020 

0.5025 

0.7066 

0.7090 

0.0021 

0.0026 

0.4999 

0.4999 

0.6971 

0.6974 

1.414 

1 .416 

0.006 

0.007 

MTBF  = 12b 
MTTR  =  1 

Gamma (2 ) 
Exponen 

0.4020 

0.4022 

0.6519 

0.6331 

0.0009 

0.0011 

0.401 1 
0.4011 

0.8273 

0.6275 

1.41 1 

1.412 

0.003 

0.004 

MTBF  =  80 
MTTR  =  2 

Gamma (2 ) 
Exponen 

1.0885 

1.0973 

1.0619 

1.0916 

0.0375 

0.0466 

1.0510 

1.0507 

0.9593 

0.9612 

2.494 

2.515 

0.08b 

0.107 

MTBF  =100 
MTTR  =  2 

Gararaa(2) 
Exponen 

0.8624 

0.8661 

0.9343 

0.9489 

0.0162 

0.0200 

0.8462 

0.8481 

0.8828 

0.8840 

2.454 

2. 4c5 

0.046 

0.057 

MTBF  =100 
MTTR  =  3 

Gamma (2 ) 
Exponen 

1 .2444 

1 .2588 

1.1474 

1.1914 

0.0598 

0.0747 

1.1846 

1.1841 

0.9977 

1.0001 

3.580 

3.o23 

0.172 

0.215 

13.  Conclusions 


Several  conclusions  are  derived  from  the  parametric  analysis: 

(a)  Restricting  the  C-teams  so  that  one  team  serves  only  a 
particular  battery  has  the  effect  of  significantly  reducing 
the  readiness  of  Battalion  units.  If  three  teams  can  serve 
any  of  3b  customer  units  on  a  dispatched,  FIFO  basis,  the 
average  number  of  nonopera tional  units  is  1.24,  versus  1.63, 
if  the  36  units  are  divided  into  3  exclusive  service  groups, 
served  by  dedicated  teams.  A  comparison  of  average  waiting 
time  in  the  service  area  is  3-38  hours,  in  the  first  case, 
versus  4.77  hours,  in  tne  second  —  a  33  %  increase. 

(b)  Little  performance  benefit  is  obtained  from  increasing  the 
average  C-team  speed  from  20  to  40  km/hr.  For  such  a  change, 
with  one  team  serving  12  customers,  the  average  number  of 
customers  down  for  this  sort  of  maintenance  cnange3  from  0.b46 
to  0.303.  Tne  average  time  spent  "in  the  maintenance  system" 
goes  from  4.77  hrs  to  4.36  hrs  witn  this  speed  doubling. 

(c)  Maintenance  performance  of  the  C-teams  suffers  a  serious 
degradation  with  a  reduction  of  the  average  speed  from  20 
km /nr  (nominal)  to  10  km /hr.  This  halving  of  the  speed 
increases  the  average  number  of  customers  in  this  service 
system  from  0.346  to  0.64b  units  (19  %  increase).  The 
average  time  spent  in  the  service  system,  i.e. , "down”,  goes 
from  4.77  hr3  to  3.71  hrs,  a  20  %  increase.  These  and 
other  parametric  effects  are  displayed  in  Tables  1,  2,  3,  and 
4  and  Figures  2  and  3. 

(d)  A  knee  occurs  at  an  MTBF  of  100  hrs  in  the  functional  relation 
of  average  customers  in  the  service  system  versus  MTBF,  for 

an  MTTR  of  3  hrs.  The  location  of  the  knee  is  somewhat 
subjective.  In  tnis  case,  a  quite  rapid  increase  occurs  in 
the  average  number  of  customers  down  for  maintenance  with  a 
small  decrease  in  the  MTBF  below  the  location  of  the  knee. 

Tne  location  of  the  knee  decreases  with  decreasing  MTTR. 
Decreasing  the  MTTR  also  makes  the  knee  more  pronounced. 

(e)  Several  alternatives  were  examined  to  describe  the  probability 
distribution  of  the  time  for  a  C-team  to  travel  to  a  customer 
ana  to  complete  maintenance  service  on  that  customer.  In 
general,  the  statistics  which  characterize  service  performance 
are  relatively  insensitive  to  the  form  of  this  distribution, 
given  a  mean  value  of  the  total  "service"  time. 
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ANNEX  A 


DERIVATION  OF  EQUATIONS 

Probability  Distribution  for  In te reus tome r  Distance 


The  cumulative  distribution  function  (c.d.f.)  for  the  distance 
a  C-team  travels  between  customers  is  derived  in  a  series  of  steps. 

We  start  by  ooserving  that  this  distance  consists  of  the  sum  of  two 
independent  random  variables,  each  of  which  is  the  projection  of 
intercustomer  distance  along  a  line.  Tnis  fact  is  a  consequence  of 
the  requirement  that  the  path  followed  by  a  C-team  consists  of  series 
of  orthogonal  segments,  each  of  whicn  is  parallel  to  one  of  the  sides 
of  the  rectangular  service  area.  One  may  equivalently  view  this  type 
of  path  as  an  x -displacement  (from  tne  customer  just  served)  followed 
by  a  y-displacement.  Let  the  dimensions  of  the  service  area  be: 
a  by  b,  so  that  0  le  x  ie  a,  0  le  y  le  b.  (Note  that  "le" 
denotes  "less  than  or  equal  to" .  Similarly,  "ge"  is  the  abbrev¬ 
iation  for  "greater  than  or  equal  to",  etc.).  Suppose  that  the 
c.d.f.  of  die  x-component  of  displacement  is  denoted  by 
F  (d  )  ana  tnat  for  the  y-component  is  F  (a  ),  wnere  the 
11  2  2 
intercustomer  distance,  d,  is 

d  =  a  +  d  ,  Gledle  a  +  b.  (1) 

1  2 

Because  the  components  of  d  are  independent,  the  c.d.f.  of  d,  F  (a), 

d 

is  the  convolution  of  F  ana  F  .  Tne  metnod  just  outlined  is 

1  2 

tne  approacn  used  to  obtain  the  distribution  of  intercustomer 
distance,  using  a  FIFO  service  discipline. 

As  a  first  step,  consider  the  distance  between  two  randomly 
selected  points  along  a  line  segment  of  unit  lengtn  (0,1),  where 
the  points  are  uniformly  distributed.  Let  the  random  location 
of  the  first  point  be  X  and  that  for  the  second  be  X  .  Note  that 

1  2 

upper  case  letters  are  used  for  realizations  of  a  random  variable. 

The  random  variable  denoting  the  distance  is  Z: 

Z  =  abs( X  -  X  ) .  (2 ; 

1  2 

Equivalently , 

Z  =  X  -  X  ,  for  X  ge  X  , 

12  12 
or 

Z  =  X  -  X  ,  for  X  ge  X  .  (j) 

2  1  2  1 
Nota tionally ,  let  tne  c.d.f.  of  z  be 

F  (z)  =  Pr(Z  le  z).  { 4 ) 

z 

Because  of  symmetry,  Pr(X  ge  X  )  is  equal  to  Pr(X  ge  X  ) . 

12  2  1 


A-1 


Therefore,  consider  only  the  first  of  tnese  two  contingencies, 
in  which  X  ge  X  . 

1  2 

Nota tionally ,  let  I  represent  the  integration  operator  over  tne 

u 

feasiole  space  of  the  variable  of  integration  u.  To  be  specific 
the  limits  of  integration  --  ul ,  u 2  —  may  be  denoted  as  arguments 
in  this  way:  I  (u1,u2).  To  abbreviate,  the  expression  (X  =  u) 
u  i 

means  that  the  random  variable  X  lies  between  u  and  u+du,  (i  =  1,2). 

i 

Then, 

F  (z)  =  I  Pr ( X  =  u)Pr(X  -  u  le  z).  (3) 

z  u  2  1 

In  performing  the  integration,  one  must  cnoose  limits  so  that 
0  le  u  le  1  and  that  X  ge  X  =  u,  as  assumed. 

1  2 

From  14)  and  (5),  remembering  that  Pr(X  =  u)  =  qu, 

i 


F  (z)  =  I  du  Pr(u  le  X  le  z+u) 
z  u  1 

F  (z)  =  I  10,1-z)  du  z  +  Pr(u  le  zj.  (b) 

z  u 

After  integrating,  evaluating  at  the  limits,  and  combining  terms, 

2 

F  (z)  =  2z  -  z  .  (7) 

z 

By  a  cnange  of  variable  witn  z  =  a  /a, 

1 

2 

F  (d  )  =  2 (d  /a)  -  Id  /a)  ,  d  le  a.  (b) 

1111  1 


Tne  probability  density  function  (p.a.f.  )  for  d 

1 

from  (b)  by  differentiation,  yielding 

f  (d  )  =  (2/a ) ( 1  -  d  /a),  d  le  a. 

11  1  1 


can  be  obtained 


(y) 


2 

Taking  the  first  two  origin  moments  of  f  —  E(d  )  and  E(d  ),  where 

1  1  1 

k  x 

E(d  )  =  I  (0,a)  u  f  (u)  du,  k  =  1,2, 

1  u  1 


A- 2 


(10) 


and 


E(d  ) 
1 


(11a) 


=  a/3 


2  2 

E(d  )  =  a  /b  .  (11b) 


The  variance  of  d  ,  V(d  ),  is  obtained  from  (11)  via  the 
1  1 

difference  between  second  origin  moment  and  mean  squared: 

2 

V  (d  )  =  a  /1 8 .  ( 1 2. 

1 

The  c.d.f.,  p.d.f.,  and  associated  moments  of  d  can  be 

2 

developed  using  the  above  results  for  d  by  merely  substituting 

1 

b  for  a  in  equations  (b)  thru  (12).  Note  that  the  random  variaole 
for  the  x-component  of  d  differs  statistically  from  that  of  tne 
y-component  only  in  its  range.  Tnus, 


2 


F  (d  )  =  2(d  /b)  -  (d  /b)  , 

2  2  2  2 

d  le  b. 

2 

(13) 

f  ( d  )  =  (2/d)  (1  -  d  /'»)  , 

2  2  2 

a  le  b . 

2 

(14) 

The  mean  and  variance  of  f  are 

2 

E(d  )  =  b/3 

2 

( 1 3a) 

and 

2 

V(d  )  =  b  / 1 8  . 

2 

(13b) 

Since  the  mean  of  the  sum  of  independent  random  variables  is 
just  the  sum  of  the  means,  from  (1),  (11a),  and  (Iba), 

E(d)  =  (a  +  b)/3.  (16) 

The  same  additive  law  applies  to  the  variance  of  the  sum  of  indep¬ 
endent  random  variables.  Tnus, 

2  2 

V (d)  =  ta  +  b  )/lb.  '17) 

Tnese  results  also  appear  on  page  4  of  this  report.  By  contrast 
to  tne  other  equations  on  page  4,  these  results  are  exact. 
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Nota tionally ,  let  tne  p.d.f.  of  d  be  denoted  by  f  (d). 

d 

As  noted  aoove,  the  probability  distribution  of  d  is  the  convol¬ 
ution  of  the  distributions  of  its  components.  This  statement 
also  applies  to  the  densities.  Both  f  (h)  and  f  (v)  are  defined 

1  2 

on  finite  domains:  (0,a)  and  (Q,b),  respectively.  Therefore, 
the  domain  of  f  (d),  witn  d  =  h+v,  is  (0,a+b). 
d 

It  will  be  seen  that  the  functional  form  changes  over  subdomains 
such  as  ( 0 ,min ( a,  b) )  and  (min ( a, b) , max( a, b) ) *  This  results  from  the 
finite  limits  of  integration  in  the  convolution  operation.  The 
specification  of  these  limits  in  the  space  of  n  and  v  is  tricky. 

A  s tra ightforward  way  of  handling  this  problem  is  to  use  Laplace 
transforms.  Tne  convolution  theorem  of  integral  transforms  applies 
here:  tne  transform  of  the  convolution  is  the  product  of  the 
transforms  of  the  functions  being  convolved.  Denoting  the  Laplace 
transform  of  f  as  f*(s),  where  i  is  1  or  2 , 
i  i 

f*(s)  =  I  (0 ,  liu( i) )  exp(-su)  f  i.u)  du, 
i  u  i 

with  liu(  1  )  =  a  and  liu(2)  =  b. 

Tnen, 

f*(s)  =  f*(s)  f*(s).  (19) 

d  1  2 

In  performing  the  inverse  transform,  the  "problem"  of  subdomains 
is  Handled  autoint ically  by  the  mechanics  of  tne  operation.  This 
attractive  aspect  of  integral  transforms  motivated  this  approach. 

After  some  algebraic  manipulation  of  the  results  from  (IB), 
the  transform  for  f  oecomes 

1 

2  2  2  2 

f*(s)  =  (2/a)/s  -  (2/a  )/s  +  (2/a  )exp(-sa)/s  .  (20) 

1 

A  comparable  expression  for  f*(s)  is  obtained  by  substituting 

2 

b  for  a  in  equation  (20).  Note  that  terms  involving  exp(-sa) 
and  exp(-sb)  appear  in  f*(s)  and  f*(s),  respectively. 

1  2 

Using  (1y),  one  notes  that  terms  involving  these  factors  and 
the  factor  expv-s(a+b))  appear  in  f*(s).  During  inversion,  terms 

d 

which  involve  a  factor  of  the  form  exp(-cs),  with  c  constant, 
are  converted  to  terms  in  the  inverse  which  have  a  unit  step 
function  as  a  factor.  Tne  unit  step  function  in  the  variable  d 
is  defined  as 

u(d  -  c)  =  0,  d  It  c 

=  1/2,  d  =  c 

=  1  f  d  gt  c.  (21 ) 

Thus,  terms  in  the  inverse  which  have  a  unit  step  factor 

contribute  only  to  the  result  for  that  portion  of  the  domain 
for  which  d  is  greater  tnan  or  equal  to  the  parameter  c. 


From  (19)  and  (20) , 


2  3 

f*(s)  =  (1/s  )A  +  (1/s  )  (B  +  B  exp(-sa)  +  B  exp(-sb))  + 
d  a  b 

4 

(1/s  )  ( C  -  C  exp(-sa)  -  C  exp(-sb)  +  C  exp( -s( a+b) ) )  , 

with 

A  =  4/(ab), 

2  2  2 
B  =  -4(a  +  D)/(ab)  ,  B  =  4/a  d  and  B  =  4/ab  , 

a  b 

2 

C  =  4/(ab)  .  (22; 

Terms  in  this  transform  involving  exp(-s(a+b))  are  needed  to 
insure  that  the  inverse  transform  (p.d.f. )  is  zero  for  d  ge  a+b. 

By  restricting  the  domain  of  the  p.d.f.  to  (0,  a+b),  these  terms 
can  be  neglected.  With  this  restriction,  the  inverse  transform 
becomes. 


2 

f  (d)  r  Ad  + (B/2)d 
d 


(B  /2)(o-a) 
a 


u(d-a)  + 


(B  A 
b 


Did 


2 

-b) 


u( d-b) 


3  3  3 

+  (C/6)(d  -  (d-a)  u(d-a)  -  (d-b)  u(d-b)),  a  le  a+b.  (23) 


From  (21)  and  (23),  it  can  be  seen  that  the  polynomial  in  d  for  the 
p.d.f.  assumes  different  forms  over  the  subdomains:  (0,  min(a,b)), 
(min(a, b; ,max(a, b) ) ,  and  (max( a, b) , a+b) .  The  c.d.f.  of  intercustomer 
distance  is  obtained  from  (23)  by  integration: 


F  (d)  =  I  (0,d)  f  (t)  dt,  a  le  a+b.  (24) 

d  t  d 

The  implementing  computer  code  for  the  result  of  ^24)  is  in  Annex  B 
on  page  B-7. 


The  c.d.f.  of  intercustomer  distance  from  (24)  was  derived 
assuming  the  coordinates  of  the  positions  of  the  two  customers 
were  statistically  independent  and  uniformly  distributed .  This 
situation  exists  under  both  FIFO  service  discipline  and  when 
less  than  2  customers  are  queued  with  SDST  discipline.  Suppose 
the  service  discipline  is  shortest  distance  (SDST).  Denote  the 
conditional  c.d.f.  of  d,  given  that  n  customers  are  waiting,  by 
G(n,d).  By  definition, 


Pr(D  gt  d,  n)  =  1  -  G(n,d).  (2b) 

But, 

Pr(D  gt  d,  n)  =  Pr(D  gt  d)Pr(D  gt  d)  ...  Pr(D  gt  d), 

1  2  n 

where  D  is  tne  random  value  of  the  n  tn  distance  from  F  (d). 
n  d 
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Thus, 


Pr(D  gt  d,  n)  =  (1  - 


F  (d)) 

Q 


and 


G(n,a)  =  1  -  (1  -  F  (d))  .  (2b; 

d 

The  unconaitional  c.d.f.  of  d  under  SDST  discipline,  G(d),  is 
obtained  by  summing,  over  n,  prouucts  of  the  form: 

PrCnumber  queued  =  n)G(d,n). 

The  probability  weights  in  the  resulting  expression  can  be  obtained 
from  the  analytic  queueing  model.  An  iterative  evaluation  of  these 
probabilities  is  required,  since  the  service  rate,  being  a  model 
input,  depends  upon  the  mean  and  variance  of  d.  The  first  two 
origin  moments  of  GCd)  can  be  obtained  by  numerical  integration 
of  the  following  equations  using  Simpson's  rule. 

E(d)  =  I  (0,d;(1  -  G(t);  dt  (27a>) 

t 

2 

E(d  )  =  I  (0, d)  t  (1  -  G(t))  dt.  (27b; 

t 

The  computer  program  which  implements  these  operations  for  any 
c.d.f.,  G(d),  is  found  in  Annex  B,  starting  on  page  B-24. 

Approximations  for  Mean  and  Variance  of  Travel  Time 


Approxima tions  are  given  on  page  4  for1  the  mean  and  variance 
of  intercustomer  travel  time.  The  derivation  of  these  results  is 
sketched  here.  The  intercustomer  time,  t,  is  the  ratio  of  the 
distance,  d,  to  the  (average)  speed,  s,  over  that  distance.  This 
(avg)  speed  is  considered  a  random  variable,  varying  from  one 
customer  to  tne  next.  Tne  mean  and  variance  of  s  are  considered 
as  given.  Since  t  is  a  function  of  d  and  s,  it  can  be  expanded 
in  a  Maclaurin  series  about  the  mean  values  of  d  and  of  s.  Tne 
mathematical  expectation  of  both  sides  of  this  equation  is  taKen. 
Cuoic-  and  higher  -terms  of  d  -  E(d)  and  of  s  -  Ets)  are  discarded. 
This  procedure  ieads  to  the  approximation  for  E(t;,  given  on  p.  4. 
The  Maclaurin  expansion  for  t  is  squared,  and  mathematical  expecta¬ 
tions  are  taxen  on  both  sides  of  this  equation.  After  discarding 
the  higher-order  terms,  one  obtains  an  approximation  for  the 
second  origin  moment  of  t.  The  approximation  for  the  variance  of 
t  is  this  second  origin  moment  minus  the  square  of  the  above 
mean.  The  neglect  of  tne  drd  and  higher  central  moments  in 
the  above  derivations  may  cause  one  to  be  concerned  about  the 
accuracy  of  these  approximations  when  the  distributions  of  d  and 
of  s  are  rather  skewed.  To  evaluate  the  error  of  approximation 
one  can  compare  tnese  results  with  exact  results  for  a  special  case. 
Tne  special  case,  uiscussed  next,  involves  distributional  assump¬ 
tions  for*  intercustomer  distance  and  speed.  Tne  c.d.f.  of  t  for 
tne  special  case  is  given  in  closed  form.  Using  this  expression, 
the  mean  and  variance  of  t  are  calculated  via  the  procedure  in  (27;. 
The  aproxiraation  has  nearly  u-digit  accuracy  for  an  example  with 
a  =  15  km  and  with  s  uniform  over  tne  interval  15  to  25  km /nr. 
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Derivation  of  the  Distribution  of  Travel  Time  for  a  Special  Case 


To  obtain  a  probability  distribution  of  intercustomer  distance, 
let  the  population  be  distributed  uniformly  along  the  x-axis,  from 
zero  to  a.  In  this  case  the  FIFO  c.d.f.  of  intercustomer  distance 
is  given  by  (b),  with  the  p.d.f.  given  by  (9).  Note  that  the 
distance  distribution  is  positively  sxewed,  and,  thus  has  a  non¬ 
zero  3rd  central  moment  --  contrary  to  the  above  implicit  assumption 
For  the  present  purpose,  denote  tne  p.d.f.  and  c.d.f.  of  the 
intercustomer  distance,  x,  as  f  (x)  and  F  (x),  respectively. 

x  x 

To  derive  the  c.d.f.  of  the  intercustomer  travel  time,  assume 
that  the  speed  nas  a  uniform  distribution  between  u  and  v. 

Thus, 

f  (s)  =  1/(v-u),  u  le  s  le  v.  (2d) 

s 

Notationally ,  let  the  c.d.f.  of  time,  t,  be  F  (t),  where 

t 

F  (t)  =  Pr(T  le  t),  for  a  random  time  T.  (29) 

t 

Observe  that  T  can  take  any  value  from  zero  to  (a/u).  For  a 
particular  realization,  T  =  X/S,  with  X  being  the  particular 
distance  and  S  Deing  the  particular  value  of  tne  intercustomer 
speed.  Then,  tne  conditions  for  T  le  t  are  : 

X  =  x  and  v  ge  S  ge  x/t,  with  0  le  x  le  a. 

Applying  these  conditions  to  (29), 


F  (t)  =  I  f  (x)  Pr(S  ge  x/t)  dx,  (30a) 

t  xx 

where 

Pr(S  ge  x/t)  =1,  x  le  ut, 

=  1  -  (x/t  -  u)/(v-u),  u  le  x/t  le  v, 

=0,  x  gt  vt.  (30b) 

Thus, 


F  (t)  =  I  (0  ,ut)  f  (x)  dx  +  I  (ut, a)  f  (x)(1  -  ( x/t-u) /( v-u) )  dx 

t  X  X  XX 

for  vt  gt  a,  (31a) 


F  (t)  =  F  tut)  +  I  (ut,vt)  f  (x)(1  -  tx/t-u) /( v-u) )  dx, 
t  x  x  x 

for  vt  le  a,  or  for  t  le  a/v.  (bib) 

Carrying  out  the  indicated  integrations  yields  the  following  results 
For  v a/u  gt  vt  gt  a,  (31a)  becomes 

2  3  2  2 

F  (t)  =  1  -  a/(3t(v-u))  +  u  t/(a(v-u))  -  (2/3)u  t  /(a  (v-u)) 
t 

+  (u/(v-u))(1  -  F  (ut)),  a/u  ge  t  gt  a/v.  (32a) 

x 

F  (t)  =  F  (vt)(1  +  u/(v-u))  -  F  (ut)  u/(v-u) 

t  X  X 

2  3  b  2 

-  (v+u)/a  t  +  2/(3a  (v-u))(v  -  u  )t  ,  t  le  a/v.  tb2b) 
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ANNEX  B 


COMPUTER  SOURCE  PROGRAMS 


The  source  programs  listed  nere  are  written  in  SIMSCRIPT  2o  for 
the  PRIME  750  minicomputer.  However,  the  code  does  not  employ  features 
unique  to  this  computer.  The  MAIN  and  subprograms  calculate  statistics 
for  the  stochastic  steady  state  of  a  maintenance  service  system.  Tnis 
system  consists  of  several  mobile  contact  teams  (C-teams;  which  are  cent¬ 
rally  dispatched.  Tne  teams  serve  a  fixed  population  of  customers  within 
a  service  area.  Tne  customers  are  assumed  to  be  uniformly  distributed 
within  a  rectangular  service  area  whose  dimensions  are  input  parameters. 
All  program  inputs  are  read  interactively,  with  prompting  messages 
displayed  at  the  terminal.  At  the  beginning  of  each  program  listing  is 
found:  (a)  a  functional  description  of  the  program,  (b)  a  list  of  the 

program  inputs  —  or  input  arguments,  for  subprograms,  (c)  a  list  of  the 
program  outputs,  and,  in  some  instances,  (d)  a  list  of  Key  enuogenous 
variables.  All  utility  functions  and  routines  are  included  in  this 
listing.  No  external  files  are  used.  Since  the  output  is  lengthy,  it 
is  necessary  to  set  up  a  COMO  file  to  display  all  of  It  and  to  obtain  a 
permanent  copy. 

The  programs  found  in  this  annex  are  located  as  follows: 


Page 

MAIN  •• CONTACT. Q 

B-2 

CDT.FIFO 

B-7 

FINITE. ME2.Q 

B-8 

SGEFA 

B-gQ 

SGEDI 

B-22 

MINRV 

B-24 

LIMSTATE 

B-2b 

B-1 


1 

2 

3 

4 

5 

b 

7 

8 

1 

2 

3 

4 

b 

D 

7 

8 

9 

10 

1 1 

12 

13 

14 

1b 

10 

17 

18 

19 

20 

21 

22 

2j 

24 

2b 

2u 

27 

28 

29 

30 

31 

32 

33 

34 

3  b 

3o 

37 

38 

39 

40 

41 

42 

43 

44 


PREAMBLE  "  CONTACT. Q 
NORMALLY  MODE  IS  REAL 

DEFINE  PDT .FIFO  AS  A  REAL  FUNCTION  WITH 

DEFINE  CDT . FIFO  AS  A  REAL  FUNCTION  WITH 

DEFINE  PDT.SDST  AS  A  REAL  FUNCTION  WITH 

DEFINE  CDT.SDST  AS  A  REAL  FUNCTION  WITH 

DEFINE  FUN. CDF  TO  MEAN  CDT. FIFO 
END  ' ' PREAMBLE 


8  ARGUMENTS 
3  ARGUMENTS 
3  ARGUMENTS 
8  ARGUMENTS 


MAIN  " CONTACT. Q 


' ' PROGRAM  SOLVES  THE  STEADY-STATE  QUEUEING  PROBLEM  ASSOCIATED  WITH  A 
' 'MOBILE  CONTACT  TEAM  'WHICH  VISITS  AND  SERVES  CUSTOMERS  IN  A  RECTANGULAR 
''SERVICE  AREA.  THE  CUSTOMERS,  REQUIRING  MAI NT  SERVICE,  ARE  ASSUMED 
''TO  BE  UNIFORMLY  DISTRIBUTED  WITHIN  A  RECTANGLE  WITH  DIMENSIONS  ARANGE 
''BY  BRANGE.  THE  POPULATION  OF  CUSTOMERS  IS  FINITE  (=  MPOP).  THE  NUMBER 
"OF  C-TEAMS  IS  EQUIVALENT  TO  THE  NUMBER  OF  SERVICE  CHANNELS  (NSERVE). 

"TWO  TYPES  OF  SERVICE  DISCIPLINE  ARE  PERMITTED:  FIFO  AND  SDST.  USING 
''FIFO,  THE  CONTACT  TEAM  TRAVELS  TO  THE  CUSTOMER  HAVING  THE  FIRST  REQUEST. 
''USING  SDST,  THE  CONTACT  TEAM  TRAVELS  TO  THE  CUSTOMER  HAVING  THE  SHORTEST 
''DISTANCE  FROM  THE  CURRENT  LOCATION.  AN  ERLANG(2)  DISTRIBUTION  IS  USED 
' 'TO  APPROXIMATE  THE  TOTAL  TIME  TO  TRAVEL  TO  AND  SERVE  A  CUSTOMER. 

t  t 

' 'INPUTS: 

' 'FLAG.DIST  _  AN  INTEGER  FLAG  TO  INDICATE  PRINTING  OF  THE  PROBABILITY 

'  '  DISTRIBUTION  OF  INTERCUSTOMER  DISTANCE. 

''FLAG. FIFO  _  AN  INTEGER  FLAG  TO  INDICATE  A  FIFO  SERVICE  DISCIPLINE 

"  U  1).  THE  ALTERNATIVE  DISCIPLINE  IS  TO  SERVE  THE  CUST- 

"  OMER  HAVING  THE  SHORTEST  DISTANCE  (SDST)  NEXT. 

' ' I PRINT  _  FLAG  TO  PRINT  QUEUEING  STATISTICS  FROM  A  SUBROUTINE. 

''ARANGE  _  THE  FRONTAL  DIMENSION  (KM)  OF  THE  RECTANGULAR  CUSTOMER 

' '  TERRITORY  SERVED  BY  A  SET  OF  CONTACT  TEAMS. 

"BRANGE _ THE  CROSS-FRONTAL  (ORTHOGONAL)  DIMENSION  (KM)  OF  THE 

' '  "  RECTANGULAR  CUSTOMER  TERRITORY. 

' 'SPEED  _  THE  AVERAGE  SPEED  (KM/HR)  OF  MOVEMENT  OF  A  CONTACT  TEAM 

' '  BETWEEN  CUSTOMERS. 

' 'MPOP  _  SIZE  OF  THE  CUSTOMER  POPULATION,  I.E.,  THE  NUMBER  OF  UNITS 

' '  SERVED  BY  CONTACT  TEAMS. 

' 'NSERVE  _  THE  NUMBER  OF  CONTACT  TEAMS  (OR  SERVICE  CHANNELS;. 

' 'HTBF  _  MEAN  TIME  (HOURS)  BETWEEN  FAILURES  OF  THE  AGGEGATE  OF  PARTS 

' '  REQUIRING  SERVICE  BY  A  CONTACT  TEAM. 

' ' M TTR  _  MEAN  TIME  TO  DIAGNOSE  AND  REPAIR  A  FAILURE  WHICH  REQUIRES 

' '  A  CONTACT-TEAM  SERVICE  REQUEST. 

I  I 

"ENDOGENOUS  VARIABLES: 

"P.NSTATE(*)  _  A  REAL  VECTOR  OF  MARKOV-STATE  PROBABILITIES  FOR  THE  SERVICE 
"  ~  SYSTEM.  P . NSTATE(N)  IS  THE  PROB  THAT  N  UNITS  ARE  IN  THE 

"  "SERVICE  SYSTEM",  I.E.,  EITHER  BEING  SERVED  OR  AWAITING 

’ '  SERVICE. 

"LAMBDA _ AVERAGE  ARRIVAL  RATE  PER  CUSTOMER  (PER  HOUR). 

"MU1 _ ^  AVERAGE  UNIT  SERVICE  RATE,  INCLUDING  TRAVEL,  (PER  HOUR). 

’ ' MLIM  _ LIMITING  NUMBER  OF  CUSTOMER  PERMITTED  IN  THE  SERVICE  SYS. 

' ' PLIM  _  PROB  ACCURACY  IN  C.D.F.,  TO  LIMIT  THE  MARKOV  STATE  SPACE. 
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4b  '  ' 

4b  ' 'OUTPUTS: 

47  ' P.NULL  _ STEADY-STATE  PROBABILITY  THAT  THE  SERVICE  SYSTEM  IS  EMPTY. 

4B  ’ ' P- SYS. FULL  _  PROBABILITY  THAT  ALL  SERVICE  CHANNELS  ARE  BUSY. 

49  P.CUST.WAIT  _  PROS  THAT  AN  "ARRIVING"  CUSTOMER,  I.E.,  A  JUST  FAILED 

50  ' '  UNIT,  MUST  QUEUE. 

''ENSYS  _  EXPECTED  NUMBER  OF  CUSTOMER  UNITS  DOWN  WITH  A  FAILURE, 

"  i  WHICH  REQUIRES  THE  USE  OF  A  CONTACT  TEAM. 

bJ  ' 'SDNSYS  _  STANDARD  DEVIATION  OF  CUSTOMER  UNITS  IN  THE  SERVICE  SYS. 

^4  "ENQ  _  EXPECTED  NUMBER  OF  CUSTOMER  UNITS  WAITING  FOR  A  SERVICE 

5b  '  '  CHANNEL. 

5o  ' ESYS. WAIT  _  EXPECTED  TIME  THAT  A  FAILED  CUSTOMER  UNIT  MUST  WAIT  UNTIL 

57  ' '  SERVICE  BY  A  CONTACT  TEAM  IS  COMPLETE. 

bo  ' ' EQ. WAIT  _  EXPECTED  WAITING  TIME  IN  THE  SERVICE  QUEUE. 

59  ’’EQ.WAIT.GQ  _  CONDITIONAL  AVERAGE  WAITING  TIME  IN  THE  SERVICE  QUEUE, 

50  ' '  GIVEN  THAT  THE  CUSTOMER  MUST  QUEUE. 

' 'SDQ.WAIT.GQ  _  CONDITIONAL  STANDARD  DEVIATION  OF  THE  WAITING  TIME  IN  THE 
u2  ' '  SERVICE  QUEUE,  GIVEN  THAT  THE  CUSTOMER  MUST  QUEUE. 

63  ”E.BUSY.SERVERS_AVERAGE  NUMBER  OF  CONTACT  TEAMS  WHICH  ARE  BUSY,  I.E., 

04  ' '  THAT  ARE  DISPATCHED  TO  A  CUSTOMER  OR  PERFORMING  SERVICE. 

bb  ”SD.BUSY.SERVERS_STANDARD  DEVIATION  OF  THE  NUMBER  OF  BUSY  SERVERS  C TEAMS ) . 

6b  ' ' 

6(  DEFINE  FLAG. DIST, FLAG. FIFO, I, IPRINT, J , K, L ,M, MLIM, MPOP, N , NSERVE  AS  INTEGER 

68  VARIABLES 

69  DEFINE  ANSWER  AS  A  TEXT  VARIABLE 

70  DEFINE  P.NSTATE  AS  A  REAL,  1 -DIMENSIONAL  ARRAY 

71  PRINT  1  LINE  THUS 

DO  YOU  WANT  FIFO  SERVICE?  (. Y  OR  N).  ALTERNATIVE  IS  SHORTEST  DISTANCE. 

73  READ  ANSWER 

74  IF  SUBSTR. F( ANSWER, 1 ,1 )  =  "Y" 

7b  LET  FLAG. FIFO=1 

76  OTHERWISE 

77  LET  FLAG.FIFOrO 

76  ALWAYS 

79  PRINT  1  LINE  THUS 

INPUT  THE  FRONTAL  DIMENSION  (KM)  OF  THE  AREA  OF  THE  CONTACT  TEAFUS) . 

8l  READ  ARANGE 

32  PRINT  1  LINE  THUS 

INPUT  THE  DEPTH  (.CROSS-FRONTAL)  DIMENSION  (KM)  OF  THE  AREA  OF  THE  TEAMvS) . 

84  READ  BRANGE 

65  PRINT  1  LINE  THUS 

INPUT  THE  AVG  SPEED  (KM/HR)  OF  THE  CONTACT  TEAM  BETWEEN  CUSTOMERS. 

87  READ  SPEED 

86  PRINT  1  LINE  THUS 

INPUT  THE  PLUS  OR  MINUS  MAX  DEVIATION  IN  SPEED  FROM  AVERAGE  (KM/HR). 

90  READ  DELTAS 

91  LET  DELTAS=MIN.F(SPEED/2.0 , DELTAS) 

92  LET  VS=DELTAS**2/;> .  0 

93  PRINT  1  LINE  THUS 

INPUT  THE  NUMBER  OF  CUSTOMERS  (FIRE  UNITS)  SERVED  BY  THE  CONTACT  TEAMS. 

9b  READ  MPOP 

96  RESERVE  P.NSTATE(*)  AS  MPOP 

97  PRINT  1  LINE  THUS 

INPUT  THE  NUMBER  OF  CONTACT  TEAMS  (SERVICE  CHANNELS)  CONSIDERED,  LE  3. 

99  READ  NSERVE 

100  LET  NSERVE=MIN.F(3 .NSERVE) 


/ 


101  PRINT  1  LINE  TH'JS 

INPUT  THE  MEAN  TIME  (HRS)  BETWEEN  PERTINENT  FAILURES  OF  A  FIRE  UNIT. 

10i  READ  MTBF 

104  PRINT  1  LINE  THUS 

INPUT  THE  AVERAGE  DIAGNOSTIC  AND  REPAIR  TIME  (HRS)  FOR  THESE  FAILURES. 
100  READ  HTTR 

107  PRINT  1  LINE  THUS 

DO  YOU  WANT  TO  PRINT  THE  PROB  DIST  OF  INTERCUSTOMER  DISTANCE?  (Y  OR  N). 
104  READ  ANSWER 

110  IF  SUBSTR. F( ANSWER, 1 , 1  )  =  "Y» 

*11  LET  FlAG.DIST=1 

112  OTHERWISE 

11b  LET  FLAG. DIST=2 

114  ALWAYS 

Ho  PRINT  1  LINE  THJS 

DO  YOU  WANT  TO  PRINT  FROM  THE  QUEUEING  SUBROUTINE?  (Y  OR  N). 

117  READ  ANSWER 

11b  IF  SUBSTR. F( ANSWER, 1,1)  =  "Y" 

119  LET  IPRINT=1 

120  OTHERWISE 

121  LET  IPRINTsO 

122  ALWAYS 

1 2 j  SKIP  3  LINES 

124  PRINT  2  LINES  THUS 

INPUT  DATA  FOR  PERFORMANCE  OF  CONTACT  TEAMS 


127  '  ' 

120  ''CALCULATE  AVG  AND  STD  DEV  OF  INTERCUSTOMER  DISTANCE  AND  TRAVEL  TIME. 
124  '  ' 

IqO  if  FLAG.FIF0=1 

131  PRINT  1  LINE  THJS 

SERVICE  DISCIPLINE  USED  IS  "FIFO". 

133  LET  ER=ARANGE/3.0+BRANGE/3.0 

134  LET  VR=(ARANGE**2+BRANGE**2)/1b.O 

1 33  OTHERWISE 

lob  PRINT  1  LINE  THUS 

SERVICE  DISCIPLINE  USED  IS  " SHORTEST-DISTANCE" . 

13b  "  LET  ER=(ARANGE+BRANGE)/3.0 

139  "  LET  VR=(ARANGE**2+BRANGE**2)*(1 .0/13.0-1 .0/23.0) 

140  CALL  MINRV  (0,2 , ARANGE+BRANGE, ARANGE,BRANGE)  YIELDING  ER, VR 

141  ALWAYS 

142  LET  ETM. TR VL=ER/ SPEED +ER# VS/ SPEED* *3  ''APPROXIMATELY 

143  LET  VTM. TRVL=V R/SPEED**2+ER**2*VS/SPEED**4+ER**2*VS**2/SPEED**B 

144  +VR*VS/SPEED**4  ''APPROXIMATELY  FROM  MACLAUREN  EXPANSION  ABOUT  AVG'S 

149  PRINT  7  LINES  WITH  AR ANGE , BRANGE, SPEED, DELTAS, MPOP , NSERVE , MTBF , MTTR 

14o  THUS 

FRONTAL  DIMENSION  vKM)  OF  SERVICE  AREA  _  **.**»» 

DEPTH  DIMENSION  (KM)  OF  SERVICE  AREA  **.**** 

RANGE  OF  SPEED  (KM/HR)  OF  CONTACT  TEAM  **.****  +  /_  «#.**»# 

NUMBER  OF  CUST  UNITS  IN  SERVICE  AREA  __  ** 

NUMBER  OF  CONTACT  TEAMS  IN  SERVICE  AREA  ** 

AVG  TIME  (HR)  BET  SERVICE  REQUESTS/CUST  ****.** 

AVERAGE  DIAGNOSTIC  AND  REPAIR  TIME  (HR)  #».**** 

134  SKIP  1  LINE 
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^bb  PRINT  7  LINES  WITH  ER, SQRT.F(VR) , SPEED, SORT. F(VS) , ETM.TRVL, 

15b  SQRT.F(VTM.TRVL),ETM.TRVL/(ETM.TRVL+MTTR) 

1 57  THUS 

AVERAGE  DISTANCE  (KM)  BETWEEN  CUSTOMERS  «».««** 

STD  DEV  DISTANCE  (KM.)  BETWEEN  CUSTOMERS  ##.**## 

AVERAGE  SPEED  (KM/HR)  BETWEEN  CUSTOMERS  *#.#*#* 

STD  DEV  SPEED  (KM/HR)  BETWEEN  CUSTOMERS  ##.#### 

AVG  TRAVEL  TIME  (HR)  BETWEEN  CUSTOMERS  **.*### 

S.D.  TRAVEL  TIME  (HR)  BETWEEN  CUSTOMERS  ««.«««* 

RATIO:  TRAVEL  TIME/( TRAVEL+REPAIR  TIMES)  «.»#*# 

^bb  SKIP  2  LINES  ~ 

166  IF  FlAG.DIST  NE  1 

167  GO  TO  LO 

168  OTHERWISE 

169  PRINT  o  LINES  THUS 

PROBABILITY  DISTRIBUTION  OF  INTERCUSTOMER  DISTANCE 


DIST  C.D.F  OF  MIN  DIST  FOR  N  CUST.  N: 

(KM)  FIFO  2  j  4 

176  LET  DR=(ARANGE+BRANGE)/20.0 

177  FOR  1=1  TO  20  DO 

178  LET  R=DR*I 

!7y  LET  CDF1=CDT.FIFO(ARANGE,BRANGE,R) 

160  LET  CDF2=1 .0-( 1 .0-CDF1 )##2 

161  LET  CDF3=1 .0-( 1 .0-CDF1 )**S 

162  LET  CDF*+=1  .0-(  1 .0-CDF1  j**4 

iy3  PRINT  1  LINE  WITH  R,CDF1 , CDF2 , CDF3 , CDF4 

184  THUS 

**.**»  *.#***  *.***x  »##**  *  **** 

18b  IF  CDF1  GE  0.9999 

187  GO  TO  KO 

168  OTHERWISE 

189  LOOP  "OVER  I 

190  'KU' PRINT  2  LINES  THUS 


193  'LO'LET  LAMBDA=1 .O/MTBF 

194  LET  MU1=1 .0/(ETrt.TRVL+MTTR) 

195  SKIP  2  LINES 

19b  PRINT  8  LINES  WITH  LAMBDA, MPOP*LAMBDA, MU1 ,NSERVE*MU1 .MPOP.NSERVE 

197  THUS 

INPUT  PARAMETERS  FOR  A  STEADY-STATE,  FINITE  wUEUEING  SYSTEM 


ARRIVAL  RATE  (LAMBDA)  _  **.*#**» 

MAX  ARRIVAL  RATE  **.***## 

SERVICE  RATE  (MU17  **.*»*»* 

MAX  SERVICE  RATE  **.***#» 

CUSTOMER  POPULATION-  ** 

SERVICE  CHANNELS  ** 


PER  HR  PER  OPERATING  CUSTOMER  UNIT 

UNITS  PER  HR 

PER  HR  PER  CUSTOMER 

UNITS  Pc,R  HR 

CUSTOMERS 

SERVERS 
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2vJu 

207 

20b 

209 

210 
211 
212 

213 

214 

21  o 
217 
21b 

219 

220 
221 
222 

223 

224 


i  i 

' 'GtT  A  LIMIT  ON  THE  SYSTEM-STATE  INDEX  (MLIM)  WHICH  SATISFIES  ACCURACY 
' ' REQUIREMENTS. 

I  I 


LET  PLIM=0 .0001 

CALL  LIMSTATE  (LAMBDA, MU1 , MPOP, NSER VE , PLIM, 0 )  YIELDING  MLIM, P. NULL, 
P.NSTATEl*) 

PRINT  1  LINE  WITH  MLIM, PLIM 
THUS 

SYSTEM  INDEX  =  **  FOR  _  #.#*»#*  ACCURACY  IN  CUM  PROB 

SKIP  H  LINES 

CALL  FINITE. MEB. Q( LAMBDA, MU  1 , MPOP, MLIM, NSERVE , IPRINT)  YIELDING  P.NULL, 
P. SYS. FULL, P. OUST . WAIT, ENSYS, SDNSYS, ENQ , ESYS. WAIT, EQ . WAIT, EQ. WAIT.GQ, 
SDU. WAIT. Gw, E. BUSY. SERVERS, SD. BUSY. SERVERS, P.NSTATEC*) 

SKIP  B  LINES 
IF  IPRINT  NE  1 

PRINT  7  LINES  WITH  P. NULL, P. SYS. FULL, ENSYS, SDNSYS, ESYS. WAIT, 

E . BUSY . SER  VERS , SD . BUSY . SERVERS 
THUS 

PROBABILITY  THAT  SERVICE  SYSTEM  IS  EMPTY  (NO  UNITS  DOWN) 

PROBABILITY  THAT  ALL  SERVICE  CHANNELS  ARE  BUSY 
AVERAGE  NUMBER  OF  CUSTOMERS  IN  THE  SERVICE  SYSTEM 
STD  DEV  NUMBER  OF  CUSTOMERS  IN  THE  SERVICE  SYSTEM 
UNCOND  AVG  WAITING  TIME  (HR)  IN  THE  SERVICE  SYSTEM 
AVERAGE  NUMBER  OF  BUSY  SERVICE  CHANNELS 
STD  DEV  NUMBER  OF  BUSY  SERVICE  CHANNELS 


ft  . 

ft  ftftft 

*. 

ftftftft 

ft  ft 

.ftftft 

** 

.  ftftft 

*  ft 

.ftftft 

ft# 

.  ftftft 

ft  ft 

.ftftft 

EBB 

ALWAYS 

EB  j 

SKIP  E  LINES 

B34 

STOP 

ebb 

END  "  CONTACT. Q 

1  FUNCTION  CDT.FIFO  (A,B,T) 

2  ' ' 

3  ''FUNCTION  EVALUATES  THE  CUM  DISTRIBUTION  FUNCTION  (C.D.F;  FOR  THE 

4  "INTERCUSTOMER  DISTANCE  OF  A  DiSPACHED  MOBILE  SERVICE  SYSTEM,  UNDER  A 

5  "FIFO  SERVICE  DISCIPLINE.  THE  CUSTOMER  AREA  (OR  TERRITORY;  IS  A  RECTANGLE 

6  ' 'WITH  DIMENSIONS  A  BY  B.  CUSTOMERS  ARE  ASSUMED  TO  BE  UNIFORMLY  DISTRI- 

7  "BUTED  IN  THIS  AREA.  INDEPENDENCE  OF  THE  X-  AND  Y-  COORDINATE  POSITIONS 

8  "OF  THE  CUSTOMERS  IS  ASSUMED.  THE  PATH  THAT  THE  SERVICE  SYSTEM  IS 

9  ' 'ASSUMED  TO  FOLLOW  IS  A  SERIES  OF  ORTHOGONAL  SEGMENTS  PARALLEL  TO  THE 

10  ' 'SIDES  OF  THE  RECTANGULAR  SERVICE  AREA. 

11  " 

12  IF  T  GT  A+B 

13  RETURN  WITH  1.0 

14  OTHERWISE 

13  LET  AB=A*B 

It)  LET  ABS=AB**2 

17  LET  CDF=2.0*T**2/AB*(1.0-T*(A+B)/AB/3.0)+T**4/ABS/G.0 

18  IF  T  GT  B 

ADD  (T-B)*»S/b.0/AB*(4.0/B-(T-B;/AB;  TO  CDF 

20  ALWAYS 

21  IF  T  GT  A 

22  ADD  (T-A)**3/6.0/AB*(4.0/A-(T-A;/AB;  TO  CDF 

23  ALWAYS 

24  RETURN  WITH  CDF 
2d  END  "CDT.FIFO 
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1 

2 

3 

4 

3 

6 

7 

a 

9 

10 

11 

12 

13 

14 

la 

Id 

17 

1b 

iy 

20 

21 

22 

2a 

24 

23 

2D 

<27 

2b 

29 

30 

31 

32 

33 

34 

3a 

36 

37 

3b 

39 

40 

41 

42 

43 

44 

43 

4o 

47 

4b 

49 

a0 

a  1 

32 

aa 

34 


ROUTINE  FINITE. ME2.Q  (LAMBDA, MU1 ,MPOP, MLIM, NSERVE , MODE)  YIELDING  P. NULL, 
P. SYS. FULL,  P. GUST. WAIT,  LNSYS,  SDNSYS,  ENQ,  E3YS.WAIT,  EO.WAIT, 

EQ. WAIT. GO,  SOU. WAIT. GO,  E. BUSY. SERVERS,  SD. BUSY . SERVERS,  P.NSTATE 


''THIS  ROUTINE  CALCULATES  THE  STEADY-STATE,  SYSTEM  STATE- PROBABILITY 
''VECTOR  FOR  A  SERVICE  SYSTEM  HAVING  EXPONENTIAL  INTERARRIVAL  TIMES  FOR 
' 'EACH  MEMBER  OF  A  FINITE  CUSTOMER  POPULATION,  AND  HAVING  NSERVE  SERVERS, 
' 'WITH  SERVICE  TIMES  DISTRIBUTED  AS  ERLANG  WITH  SHAPE  PARAMETER  2. 


’ ' INPUT: 

' 'LAMBDA 
'  'MU1 
'  ' MPOP 
' 'MLIM 
' 'NSERVE 
' 'MODE 


ARRIVAL  RATE  PER  INDIVIDUAL  IN  THE  POPULATION 
SERVICE  RATE  OR  RECIPROCAL  OF  MEAN  SERVICE  TIME 
CUSTOMER  POPULATION  SIZE 

MAX  ALLOWED  CUSTOMERS  IN  THE  SERVICE  SYSTEM 
NUMBER  OF  SERVERS  (N  LE  5) 

INTEGER  FLAG  FOR  PRINTING  FROM  ROUTINE  (FOR  MODE-1 ) 


' 'OUTPUT: 


' 'P.NSTATE 


' ' P. NULL 
' ' P. SYS. FULL 
' ' P. GUST. WAIT 
'  'ENSYS 
' ' SDNSYS 
'  '  ENQ 

' ' ESYS. WAIT 
' 'EO.WAIT 
' ' EG . WAIT. GQ 


VECTOR  OF  PROBABILITY  ELEMENTS  IN  WHICH  THE  N  TH 
ELEMENT  IS  THE  PROBABILITY  THAT  N  INDIVIDUALS  ARE 
IN  THE  SERVICE  SYSTEM  IN  THE  STEADY  STATE. 
PROBABILITY  THAT  THE  SYSTEM  IS  EMPTY 
PROBABILITY  THAT  ALL  SERVICE  CHANNELS  ARE  FULL 
PROBABILITY  THAT  AN  ARRIVING  CUSTOMER  MUST  QUEUE 
EXPECTED  NUMBER  OF  CUSTOMERS  IN  THE  SYSTEM 
STD  DEVIATION  OF  THE  NUMBER  IN  THE  SYSTEM 
EXPECTED  NUMBER  OF  CUSTOMERS  IN  THE  UUEUE 
EXPECTED  WAIT  IN  THE  SERVICE  SYSTEM 
EXPECTED  WAIT  IN  THE  SERVICE  QUEUE 
EXPECTED  QUEUE  WAITING  TIME,  GIVEN  THAT  A 
CUSTOMER  MUST  QUEUE. 


' ' SDQ. WAIT.GQ  STANDARD  DEVIATION  OF  QUEUE  WAITING  TIME,  GIVEN 

' '  THAT  A  CUSTOMER  MUST  QUEUE. 


'' E. BUSY  .SERVERS  EXPECTED  VALUi.  JF  BUSY  SERVERS 

' 'SD. BUSY. SERVERS  STANDARD  DEVIATION  OF  THE  NUMBER  OF  BUSY  SERVERS 

' 'OPTIONAL  PRINTING  OF  THE  UUEUE  WAITING  TIME  DISTRIBUTION  IS  PROVIDED. 


' ' ENDOGENOUS 

’  'AM 
'  'BV 

!  I 
t  t 
f  I 

1  *gv 

i  t 
t  t 
i  i 
t  i 

•  1 IPVTC* ) 

i  i 

'  ' DETt  * ) 

t  I 
t  I 


VARIABLES: 

KOLMOGOROV  ST ATE- TRANSITION  MATRIX 

VECTOR  OF  PROBABILITIES  IN  WHICH  THE  N  TH 

ELEMENT  IS  THE  PROBABILITY  OF  FINDING  THE 

SYSTEM  IN  THE  N  TH  STATE,  WHERE 

N=2* (NO  CUSTOMERS- 1 )+ST AGE  OF  SERVICE,  NO  GT  0. 

VECTOR  OF  PROBABILITIES  IN  WHICH  THE  N  TH 
ELEMENT  IS  THE  PROBABILITY  THAT  AN  ARRIVING  CUSTOMER 
FINDS  THE  SYSTEM  IN  THE  N  TH  STATE.  THIS  VECTOR 
IS  USED  TO  CALCULATE  THE  DISTRIBUTION  OF 
WAIT  IN  QUEUE. 

INTEGER  VECTOR  USED  IN  FACTORING  THE  STATE-TRANSITION 
MATRIX. 

TWO- ELEMENT  VECTOR  USED  TO  STORE  THE  DETERMINANT  OF 
THE  STATE-TRANSITION  MATRIX. 


DEFINE  I, INFO, J, MLIM, MPOP, TWOM, MAX, MAXM, MODE, N, NSERVE  AS  INTEGER  VARIABLES 
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55  DEFINE  IPVT  AS  AN  INTEGER,  1 -DIMENSIONAL  ARRAY 

5b  DEFINE  BV , DET , PV , QV, P. NSTATE, AND  Q.NSTATE  AS  REAL, 1 -DIMENSIONAL  ARRAYS 

57  DEFINE  AM  AS  A  REAL,  2-DIMENSIONAL  ARRAY 

58  IF  MLIM  LE  0  OR  MLIM  GT  MPOP 

59  PRINT  1  LINE  WITH  MPOP, MLIM 

60  THUS 

ERROR  IN  INPUT  TO  FINITE. ME2.Q.  MPOP  =  ***  MLIM  =  ***. 

62  STOP 

63  OTHERWISE 

64  LET  TW0M=2*MLIM 

65  RESERVE  DET(*)  AS  2 

ob  LET  MU=2.0#MU1  ''SERVICE  RATE  FOR  EACH  OF  THE  2  STAGES 

67  LET  ML=REAL.F(MPOP)#LAMBDA  ' 'MAX  ARRIVAL  RATE 

68  RESERVE  P.NSTATEC* )  AS  MPOP 

69  RESERVE  Q. NSTATE(# )  AS  MPOP 

70  IF  NSERVE  LE  0 

71  RETURN 

72  OTHERWISE 

73  IF  NSERVE  GT  1 

74  GO  TO  LO 

75  OTHERWISE  "NSERVE  =  1 

77  RESERVE  PV(*)  AS  TWOM 

78  RESERVE  QV(*)  AS  TWOM 

79  RESERVE  AM(*,*)  AS  TWOM  BY  TWOM 

80  RESERVE  IPVT(*)  AS  TWOM 

81  FOR  1=1  TO  TWOM  DO 

82  LET  QV(I)=0.0 

83  FOR  J=1  TO  TWOM  DO 

84  LET  AM(I,J)=0.0 

85  LOOP  "OVER  J 

86  LOOP  "OVER  I  TO  INITIALIZE 

87  '  ' 

88  "FILL  NON-ZERO  ELEMENTS  OF  THE  REDUCED  STATE-TRANSITION  MATRIX  (AM) 

89  "AND  THE  CONSTANT  VECTOR  (BV). 

90  '  ' 

91  FOR  J=1  TO  TWOM  DO 

92  LET  AM( 1 , J)=-ML 

93  LOOP  "OVER  COLUMNS 

94  SUBTRACT  MU+REAL.F(MP0P-1 )*LAMBDA  FROM  AM( 1,1) 

95  ADD  MU  TO  AM( 1 ,4) 

96  LET  AM(2 , 1 ) =MU 

97  LET  AM(2 ,2)=-MU-ML+LAMBDA 

98  FOR  N=2  TO  MLIM-1  DO 

99  LET  I=2#N-1 

100  LET  AM(I,I-2)=(MP0P-N+1 ) "LAMBDA 

101  LET  AM( I , I) =-MU-( MPOP-N) *LAMBDA 

102  LET  AM( 1,1+3) =MU 

103  LET  AM( 1+1 ,1-1 ) =AM( I , 1-2) 

104  LET  AM(I+1 ,I)=MU 

105  LET  AM(I+1 ,1+1 )=AM(I,I) 

106  LOOP  ' 'OVER  N 

107  "FINALLY, 

108  LET  AM(TW0M-1 ,TW0M-3)=LAMBDA*(MP0P-MLIM+1 ) 

109  LET  AM( TWOM- 1 , TWOM- 1 )=-MU 

110  LET  AM(TW0M,TW0M-2)=LAMBDA*(MP0P-MLIM+1 ) 

111  LET  AM( TWOM, TWOM-1 )=MU 

112  LET  AM( TWOM, TWOM) =-MU 
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113  " 

114  "OBTAIN  THE  INVERSE  OF  AMC*,*)  IN  PLACE. 

115  " 

11b  '  ' 

117  "FACTOR  AM(*,#;  FIRST  TO  OBTAIN  AN  INVERSE. 

118  " 

119  CALL  SGEFA  ( AM( *, * ) , IPVT(* ) , INFO) 

120  IF  INFO  NE  0 

1^1  PRINT  1  LINE  WITH  INFO  THUS 

TROUBLE  FACTORING  AM.  INFO  =  *. 

123  STOP 

124  OTHERWISE 

125  CALL  SGEDI  (AM(* ,*) , IPVT(») , DET(») , 1 1 ) 

1  26  " 

127  "OBTAIN  THE  SOLUTION  VECTOR  OF  THE  MATRIX  EQUATION. 

128  " 


129 

130 

131 

132 

133 

134 
1  3b 
136 

137 

138 

139 

140 

141 

142 

143 

144 

145 
14b 

147 

148 

149 


FOR  1  =  1  TO  TWOM,  LET  PV ( I)  =-ML*AMU  ,  1  ) 

I  I 

"SOLVE  FOR  THE  EMPTY-SI  STEM  PROBABILITY  (P.NULL)  AND  CALCULATE 

"THE  STATE-PROBABILITY  VECTOR  (.P.NSTATE). 

1  » 

LET  ENSYSsU.U 
LET  VNSYS=0.0 
LET  ENQ=0 .0 
LET  ESQ =0.0 
LET  SUM =0.0 
FOR  N=1  TO  MLIM  DO 
LET  I=2*N-1 

LET  P.NSTATE(N)=PV(I)+PV(I+1) 

ADD  P.NSTATE(N)  TO  SUM 
ADD  N*P.NSTATE(N)  TO  ENSYS 
ADD  N**2*P.NSTATEtN)  TO  VNSYS 
ADD  (N-1 )*P.NSTATE(N)  TO  ENQ 
ADD  (N-1 )**2*P.NSTATE(N)  TO  ESQ 
LOOP  "OVER  NON-ZERO  SYSTEM  STATES 

I  I 

"CHECK  VALIDITY  OF  PROBABILITY  SUM. 


150  " 

151  IF  SUM  LT  0.0  OR  SUM  GT  1.0 

1 52  PRINT  1  LINE  WITH  SUM  THUS 

ERROR  IN  ROUTINE  FINITE. ME2. Q.  PARTIAL  SUM  OF  STATE  PROBS  =  **.**#. 
154  STOP 

lb 5  OTHERWISE 


15b  LET  P.NULL=1 -O-SUM 

157  LET  P.SYS.FULL=SUM  "FOR  A  SINGLE  SERVER 

158  " 

159  "CALCULATE  THE  VARIANCE  AND  STD  DEV  OF  THE  NUMBER  IN  THE  SYSTEM. 
1  bO  '  ' 

1b 1  LET  VNSYS= VNSYS- ENSYS* *2 

162  LET  SDNSYS=SQRT.F( VNSYS) 

Ibi  LET  VNQ=ESQ-ENQ**2 

1b4  LET  SDNQ=SQRT.F(VNQ) 

165  '  ' 

1  bo  "CALCULATE  THE  MEAN  AND  STANDARD  DEVIATION  OF  BUSY  SERVERS. 

1o7  '  ' 

168  LET  E. BUSY. SERVERS=1 .0-P. NULL  "FOR  A  SINGLE  SERVER 

169  LET  VAR. BUSY. SERVERS=1 .U-P. NULL-E. BUSY. SERVERS**2 

170  LET  SD. BUSY. SERVERS=SQRT.F( VAR. BUSY. SERVERS) 
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171 

17* 

173 

174 
17b 
176 

177 

178 

179 

180 
181 
182 

183 

184 
18b 
186 

187 

188 

189 

190 

191 

192 

193 

194 
19b 

19b 

197 

198 

199 

200 
201 
202 

203 

204 

205 

206 

207 

208 

209 

210 


''CALCULATE  MEAN  WAITING  TIMES  USING  LITTlE'S  FORMULA. 

I  I 

LET  ESY S .  WAI T=EN  SY S/ LAMBDA/ (  REAL.F(MPOP)-ENSYS) 

LET  EQ.WAIT=ESYS. WAIT-1. 0/MU1 
IF  MPOP=1 
RETURN 

OTHERWISE  "CALCULATE  AND  PRINT  WAITING-TIME  DISTRIBUTIONS 

I  1 

"CALCULATE  THE  PROBABILITY  THAT  AN  ARRIVAL  FINDS  THE  SYSTEM 
"IN  THE  N  TH  STATE  (QV(N)). 

1  I 

LET  NORM. CONST =MPOP*P. NULL 
FOR  N=1  TO  MLIM-1  DO 

ADD  ( MPOP-N)*P.NSTATE(N)  TO  NORM. CONST 
LOOP  ’ ' TO  CALCULATE  THE  NORMALIZATION  CONSTANT 
LET  Q.NULL=MPOP*P. NULL/NORM. CONST 
LET  P . CUST . WAIT=1 . 0-Q. NULL 
LET  EQ. WAIT . GQ=0 .0 
LET  VQ. WAIT.GQ=0.0 
FOR  N= 1  TO  MLIM-1  DO 
LET  I=**N-1 

LET  QV(I)=(MPOP-N)*PV (I)/NORM. CONST 
LET  QV(I+1)=(MPOP-N;*PV(I+1)/NORM. CONST 
LET  Q. NSTATE (N) =QV (l)+QV(I+1 ) 

ADD  ( 1+1 )*QV(I)+I*QV(I+1 )  TO  EQ. WAIT.GQ 
ADD  (1+1 )*(I+2)*QV(I)+I*(I+1 )*QV(I+1 )  TO  VQ. WAIT.GQ 
LOOP  "OVER  THE  NUMBER  OF  ARRIVING  CUSTOMERS 
LET  EQ. WAIT. CHECK =EQ. WAIT.GQ. /MU 
LET  ESYS. WAIT. CHECK=EQ. WAIT. CHECK  +1.0/MU1 
LET  EQ. WAIT. GQ=EQ. WAIT. GQ/MU/d  .0-Q. NULL) 

LET  VQ. WAIT.GQ=VQ. WAIT.GQ /MU/MU/ ( 1 . 0-Q. NULL;-EQ. WAIT.GQ*** 

LET  SDQ . WAIT . GQ  =SQRT . F ( VQ . WAIT . GQ ) 

IF  MODE  NE  1 
RETURN 
OTHERWISE 
'Ll 'SKIP  4  LINES 

PRINT  21  LINES  WITH  NSERVE,P. SYS. FULL, ENSYS,SDNSYS,ENQ,SDNU,E. BUSY. SERVER 
SD. BUSY. SERVERS, ESYS. WAIT, EQ. WAIT, EQ. WAIT. GQ,P.NULL,P.NULL,Q.NULL,Q. NULL 
THUS 

STATE  OF  THE  SYSTEM  IN  STEADY  STATE  WITH  GAMMA(2)  SERVICE  FOR  **  SERVERS 


PROBABILITY:  SYSTEM  IS  FULL  *.***» 

EXPECTED  NUMBER  IN  THE  SYSTEM  *#.**#* 

STD  DEV  NUMBER  IN  THE  SYSTEM  **.**** 

EXPECTED  NUMBER  IN  THE  QUEUE  **.**** 

STD  DEV  NUMBER  IN  THE  QUEUE  **.**** 

EXPECTED  NUMBER  BUSY  SERVERS  **.*#*» 

STD  DEV  NUMBER  BUSY  SERVERS  **.*»»* 

MEAN  WAITING  TIME  (HRS)  IN  SYSTEM  ***.*** 

MEAN  WAITING  TIME  (HRS)  IN  QUEUE  ***.*** 

COND'AL  MEAN  TIME  (HRS)  IN  QUEUE  ***.*** 


PROBABILITY  DISTRIBUTIONS  FOR  NUMBER  OF  CUSTOMERS  IN  THE  SYSTEM 


NO  IN 
SYSTEM 


UNCONDITIONAL 

PROB  CUML 

DENS  PROB 


GIVEN  ARRIVAL 
PROB  CUML 
DENS  PROB 


232 
233 
2j  4 
235 
2  36 

237 

238 

240 

241 

242 

243 

244 


0  #.##*#  *.****  #.***#  »,«•«« 

LET  CUM. DIST=P. NULL 
LET  CUM. DIST . COND-Q .NULL 
FOR  N=1  TO  MLIM  DO 

ADD  P.NSTATE(N)  TO  CUM. DIST 
ADD  Q.NSTATE(N)  TO  CUM. DIST . COND 

PRINT  1  LINE  WITH  N,  P.NSTATE(N),  CUM. DIST,  Q.NSTATE(N) , 


THUS 

*  *  #.**## 

*,***» 

£.**** 

*  ,  *  *** 

IF  CUM. DIST  GE 

0.9999 

GO  TO  K1 


OTHERWISE 

LOOP  "OVER  THE  SYSTEM  STATES 
' K1 ' PRINT  2  LINES  THUS 


CUM. DIST. COND 


247  IF  NSERVE  GT  1 

24 8  RETURN 

249  OTHERWISE 

250  LET  DELT=MAX. F( 1 .0 , TRUNC. F (EQ. WAIT.GQ/ 10 .0 ) ) 

231  LET  MAX=o0 

252  SKIP  2  LINES 

2^3  '  ' 

254  "PRINT  HEADINGS  FOR  THE  QUEUE  WAITING  TIME  DISTRIBUTION. 

253  '  ' 

250  PRINT  8  LINES  WITH  MPOP  AND  U. NULL  THUS 

CUM  AND  CONDITIONAL  CUM  PROB  DISTRIBUTIONS  OF  WAITING  TIME  IN  QUEUE 
EXPONENTIAL  INTERARRIVALS  FOR  EACH  OF  **  CUSTOMERS. ERLANG(2 )  SERVICE. 


TIME  CUML 
(HRS)  PROB 

COND 

PROB 

o  *,*### 

0 

2o5 

t  t 

206 

' 'START  TIME  LOOP. 

267 

1  t 

208 

FOR  1=1 

TO  MAX  "TIME  STEPS"  DO 

269 

LET 

TM=I*DELT 

270 

LET 

ARG=TM*MU 

271 

LET 

CUM. DIST=1. O-EXP. F(-ARG)*(QV(1 )*(1 

272 

LET 

SUMN=0 .0 

273 

FOR 

N=2  TO  MLIM-1  DO 

274 

LET  SUMEJ  =  1.0 

275 

LET  ARGJ=1.0 

276 

LET  JFACT0RIAL=1 .0 

2  77 

FOR  J=1  TO  2*(N-1)  DO 

278 

LET  ARGJ=ARGJ*ARG 

279 

LET  JFACTORIAL^J  FACTORIALS 

280 

ADD  ARGJ/ JFACTORIAL  TO  SUMEJ 

281 

LOOP  ' 'OVER  J 

282 

LET  ARGJ=ARGJ*ARG 

285 

LET  JFACT0RIAL=JFACT0RIAL*(2*N-1 ) 

284 

LET  SUMOJ=SUMEJ +ARGJ/ JFACTORIAL 

285 

ADD  QV(2»N-1 )*SUM0J+QV(2*N)»SUMEJ  ' 

20b 

LOOP 

’  ' 'OVER  N 

28  7 

LET 

CUM. DIST =CUM. DIST-EXP.F(-ARG) *SUMN 

TO  SUMN 
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<£88 

289 


LET  CON. D 1ST =( CUM. DIST-Q. NULL)/(1 .O-U. NULL) 

PRINT  1  LINE  WITH  TM,  CUM.DIST,  AND  CON.DIST  THUS 
##*.#  *.**#*  *.#*## 

2 yl  IF  CUM.DIST  GE  0.9999 

292  GO  TO  K2 

298  OTHERWISE 

292  LOOP  ' 'OVER  I  TIME  STEPS 

299  'K2' PRINT  10  LINES  WITH  EQ.WAIT.GQ,  SDQ. WAIT.GQ,  EQ. WAIT. CHECK, 

29o  ESYS. WAIT. CHECK,  Q.NULL,  E. BUSY. SERVERS,  SD. BUSY. SERVERS 

297  THUS 


MEAN  WAITING  TIME  IN  QUEUE,  GIVEN  A  WAIT  *#*.*## 
STD  DEV  WAITING  TIME  IN  QUEUE,  GIVEN  A  WAIT  »*».##» 
MEAN  WAITING  TIME  IN  QUEUE  ( UNCONDITIONAL)  »*#.### 
MEAN  WAITING  TIME  IN  SYS  (UNCONDITIONAL)  *»».*»# 
PROBABILITY  THAT  AN  ARRIVAL  FINDS  SYS  EMPTY  ».#*»* 
MEAN  NUMBER  OF  BUSY  SERVERS  #*.##*# 
STD  DEV  NUMBER  OF  BUSY  SERVERS  ##.#### 


308  RETURN 

309  'LO'IF  NSERVE  GT  2 

310  GO  TO  L3 

311  OTHERWISE 

312  '  ' 

313  ''RESERVE  ARRAYS  FOR  THE  CASE:  NSERVE=2 

314  '  ' 

319  LET  MAXM=3*MLIM-1 

31o  "  RESERVE  BV(*;  AS  MAXM 

317  RESERVE  PV(*)  AS  MAXM 

318  RESERVE  QV(*)  AS  MAXM 

319  RESERVE  AM ( * , * )  AS  MAXM  BY  MAXM 

320  RESERVE  IPVT(»)  AS  MAXM 

321  '  ' 

322  ''FILL  THE  REDUCED  STATE-TRANSITION  MATRIX  (AM). 

323  '  ' 

324  FOR  1=1  TO  MAXM  DO 

325  FOR  J=1  TO  MAXM  DO 

32o  LET  AM( I , J ) =0.0 

327  LOOP  ' 'OVER  J 

328  LOOP  ''OVER  I 

329  FOR  J=1  TO  MAXM  DO 

330  LET  AM(1,J)=ML 

331  LOOP  "OVER  J 

332  LET  AM(1  ,2>=AMU  ,29+MU 

333  LET  AM(2 , 1 )=MU 

334  LET  AM(2,2)=-(MP0P-1 )*LAMBDA-MU 

339  LET  AM(2,5)=2.0#MU 

33t>  LET  AM(3,1  )  =  (MP0P-1  J*LAMBDA 

337  LET  AM(3,3)=-(MP0P-2;*LAMBDA-2.0»MU 

338  LET  AM( 3,7) =MU 

339  LET  AM(4,2)=(MP0P-1 1*LAMBDA 

340  LET  AM(4 ,3)=2. 0*MU 

341  LET  AM( 4,4) =AM( 3,3) 

342  LET  AM(4,8;=2.0»MU 

343  LET  AM(5,4;=MU 

344  LET  AM(  5,5)=AMt,  4,4) 
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345 

34b 

347 

348 

349 

350 

351 

352 
3b3 

354 

355 
350 
357 

3b8 

359 

360 
3o1 
362 
363 

804 

365 

3ob 

3b7 

3b8 

369 

3/0 

371 

372 

373 

374 

375 

376 

377 
376 
379 


3o1 

362 

363 

364 

365 
386 
367 

j88 

369 

390 

3^1 

392 

393 

394 

395 
39b 
397 
39b 

399 


FOR  N=3  TO  MLIM-1  DO 
LET  I=3*CN-1J 

LET  AM(I,8*N-6J=CMP0P-N+1 ) “LAMBDA 
LET  AMC I ,  I )  =-(  ( MPOP-N)  *L AMBDA+2 . 0 *MU) 

LET  AM( I,3*N+1 )=MU 
ADD  1  TO  I 

LET  AMC I,3*N-5)  =  CMPOP-N+1 1*LAMBDA 
LET  AM( I ,3*N-3)=2. U*MU 
LET  AM( I , I )  =AMl 1-1 ,1-1 ) 

LET  AM(I,3*N+2)=2.0*MU 
ADD  1  TO  I 

LET  AM( I,3*N-4jr tMPOP-N+1 1*LAMBDA 
LET  AM( I , 3*N-2) =MU 

LET  AMlI,I)=-C  wMP0P-N)*LAMBDA+2.0*MU1 
LOOP  ’’OVER  NUMBER  OF  CUSTOMERS  IN  THE  SYSTEM 
LET  I=3*CMLIM-1 ) 

LET  AMCl,3*MLIM-b;=LAMBDA*(MP0P-MLIM+1 ) 
LET  AM(I,I)=-2.0*MU 
ADD  1  TO  I 

LET  AMCI,3*MLIM-5>=LAMBDA*CMP0P-MLIM+1 ) 
LET  AM( I , 3*MLIM-3 ) =2 .  0*MU 
LET  AMCl,i;=-2.0*MU 
ADD  1  TO  I 

LET  AMC  1 , 3 *MLIM-4  J  =  LAMBDA* »,MP0P-MLIM+1 ) 
LET  AM( I , 3*MLIM-2 ) =MU 
LET  AMCI,I)=-2.0*MU 

I  I 

'  'OBTAIN  THE  INVERSE  OF  AMC  * , *  1  IN  PLACE. 


'  'FACTOR  AMC  * , * ) . 

I  I 

CALL  SGEFA  C AM( * , * ) , IPVTC * ) , INFO) 

IF  INFO  NE  0 

PRINT  1  LINE  WITH  INFO  THUS 
TROUBLE  FACTORING  AM.  INFO  =  *. 

STOP 

OTHERWISE 

CALL  SGED1  lAM( *,*), IPVTC*) , DETl*) , 1 1 ) 

I  I 

"OBTAIN  THE  SOLUTION  VECTOR  OF  THE  MATRIX  EQUATION. 

I  I 

FOR  1=1  TO  MAXM,  LET  PV  ( I )  =ML*AM(  1 , 1 ) 

I  I 

"SOLVE  FOR  THE  EMPTY-SYSTEM  PROBABILITY  CP. NULL;  AND  CALCULATE  THE 
"CUSTOMER  STATE  PROBABILITY  VECTOR  (P.NSTATE). 

?  I 

LET  ENSYS=PVC1 )+PVC2) 

LET  P.NSTATEC 1 ) =ENSYS 
LET  VNSYS=ENSYS 
LET  ENQ=0 . 0 
LET  ESQ =0.0 

LET  SUM=ENSYS  ' 

FOR  N=2  TO  MLIM  DO 
LET  1  =  3*  C  N-1 ) 
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400  LET  P.NSTATE(N)=PV(I)+PVvI+1 )+PV(I+2) 

401  ADD  P.NSTATE(N)  TO  SUM 

402  ADD  N*P.NSTAT£(N)  TO  ENSYS 

403  ADD  iN-2)*P.N3TATE(N)  TO  ENQ 

404  ADD  (N-2)**2*P.NSTATE(N)  TO  ESQ 

40s  ADD  N**2*P.NSTAT£(N)  TO  VNSYS 

400  LOOP  ’’OVER  THE  OTHER  NON-ZERO  SYSTEM  STATES 
407  '  ' 

406  ’’CHECK  VALIDITY  OF  THE  PROBABILITY  SUM. 

409  '  ' 

410  IF  SUM  LT  0.0  OR  SUM  GT  1.0 

411  PRINT  1  LINE  WITH  SUM  THUS 

ERROR  IN  ROUTINE  FINITE. MEE. Q.  PARTIAL  SUM  OF  STATE  PRObS  =  **.*»*» 

413  STOP 

414  OTHERWISE 

416  LET  P.NULL=1 .0-SUM 

41b  LET  P.SYS.FULL=1 .0-P.NJLL-P.NSTATE( 1 )  ’’FOR  2  SERVERS 

417  " 

416  ' 'CALCULATE  THE  VARIANCE  AND  STANDARD  DEV  OF  THE  NUMBER  IN  THE  SYSTEM. 

419” 

4E0  LET  VNSYS=VNSYS-ENSYS**E 

421  LET  SDNSYS=SQRT.F( VNSYS) 

422  LET  VNQ =ESQ -ENQ **2 

423  LET  SD  NQ  =SQ  RT . F ( VNQ ) 

424  '  ' 

426  ' 'CALCULATE  THE  MEAN  AND  STANDARD  DEVIATION  OF  BUSY  SERVERS. 

42b  '  ' 

427  LET  E.BUSY. SERV£RS=P.NSTATE( 1 )+2. 0* ( 1 . 0-P. NULL-P.NSTATE( . ) ) 

428  LET  VAR  .BUSY.  SERVE  RS=P.  N  ST  ATE  ( 1  )  +  4.0*U  .  0-P.  NULL-  P.N  STATU  1  ) ) 

429  -E . BUSY. SER VERS**2 

430  LET  SD. BUSY. SEHVERS=SQRT.F( VAR. BUSY. SERVERS) 

431  '  ' 

432  ' 'CALCULATE  MEAN  WAITING  TIME  USING  LITTLE'S  FORMULA. 

433  ' ' 

434  LET  ESYS. WAIT=ENSYS/LAMBDA/( REAL. F(MPOP) -ENSYS) 

435  LET  EQ. WAIT  =ENQ/LAMBDA/( REAL.F(MPOP) -ENSYS) 

43b  IF  MPOP  LE  2 

437  RETURN 

438  OTHERWISE  ' 'CALCULATE  AND  PRINT  CONDITIONAL  PROB  DISTRIBUTIONS 

439  ” 

440  ' 'CALCULATE  THE  PROBABILITY  THAT  AN  ARRIVAL  FINDS  THE  SYSTEM  IN  STATE  N. 

441  '  ' 

442  LET  NORM. CONST =MPOP*P.NULL 

443  FOR  N= 1  TO  MLIM-1  DO 

444  ADD  (MPOP-N)*P.NSTATE(N)  TO  NORM. CONST 

445  LOOP  ''TO  CALCULATE  THE  NORMALIZATION  CONSTANT 

44b  LET  Q.NULL=MPOP*P. NULL/NORM. CONST 

447  LET  QV (1 ) =(MP0P-1 )#PV ( 1 ) /NORM. CONST 

448  LET  QV(2)  =  (MP0P-1 )*PV (2) /NORM.  CONST 

449  LET  Q.MSTATE (1 ) =QV (1 ) +QV (2) 

450  FOR  N=2  TO  MLIM-1  DO 

451  LET  I=3*(N-1) 

452  LET  QV(I)= (MPOP-N)*PV(I) /NORM. CONST 

453  LET  QV (1+1 )=(MP0P-N)*PV(I+1 )/NORM. CONST 

454  LET  QVCI+2)=( MPOP-N; #PV ( 1+2 ) /NORM. CONST 

455  LET  Q. NSTATE (N ) =QV (I)+QV(I+1) +QV ( 1+2 ) 

458  LOOP  "OVER  NUMBER  OF  ARRIVING  CUSTOMERS 
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4t>7  LET  P.  GUST.  WAIT=  1  . 0-g.NULL-Q.NSTATE  ( 1 ) 

458  ' L2' LET  EO. WAIT. GQ=EQ. WAIT/P. CUST. WAIT 

43y  LET  SDQ. WAIT. GQ=EQ. WAIT. GQ  ''APPROXIMATELY 

460  IF  MODE  NE  1  RETURN 

401  OTHERWISE 

402  GO  TO  LI 

463  ' L3 ' IF  MPOP  LE  b 

404  PRINT  1  LINE  WITH  MPOP  AND  NSERVE  THUS 

ERROR  IN  FINITE. MEE. Q.  MPOP  =  ** .  NSERVE  =  ** . 
4bo  STOP 

467  OTHERWISE 

Lob  '  ' 

4oy  ' 'RESERVE  ARRAYS  FOR  THE  CASE:  NSERVE  =  3. 

470  '  ' 

471  LET  MAXM=4*MLIM-3 

473  RESERVE  PV(*)  AS  MAXM 

4/4  RESERVE  gV^*;  AS  MAXM 

47b  RESERVE  AM(*,*)  AS  MAXM  BY  MAXM 

47 o  RESERVE  IPVT(*)  AS  MAXM 

477  '  ' 

47 o  ''FILL  THE  REDUCED  STATE  TRANSITION  MATRIX  l AM) 

479  " 

48 0  FOR  1=1  TO  MAXM  DO 

482  FOR  J=1  TO  MAXM  DO 

46b  LET  AM( I, J)=0.0 

464  LOOP  ' 'OVER  J 
4bb  LOOP  ''OVER  I 

4bo  FOR  J=1  TO  MAXM  DO 
407  LET  AMU,J;=ML 

488  LOOP  ''OVER  J 

489  LET  AM( 1 ,2)=AM( 1 ,2 j+MU 

4yO  LET  AM  ( 2 , 1 )=MU 

491  LET  AM(2  ,2)=-{.  (MPOP-1  )#LAMBDA+MU) 

492  LET  AM (2  ,  b)=2.0*M'J 

493  LET  AM( 3, 1 )= (MPOP-1 )*LAMBDA 

494  LET  AM(.3,3)=-(  t.MP0P-2;*LAMBDA+2. 0*MU) 

495  LET  AM(3,7)=MU 

49b  LET  ARK  4 , 2 ;  =  ( MPOP- 1 ) *  LAMBDA 

497  LET  AM(4 ,3)=2.0*MU 

4yb  LET  AM(4,4)=AM(3,3J 

499  LET  AM(4,8)=2.0*MU 

bOO  LET  AMlb,4)=MU 

bOI  LET  AM(b,b;=AM(4,A) 

bOE  LET  AM(b,9)=3.0*MU 

bOj  LET  AMC  o , 3 )  =  ( MPOP-2 ) *LAMBDA 

b04  LET  AM(6,b)  =  -UMP0P-b)#LAMBDA+b.0*MU) 

bOb  LET  AM(6 , 1 1 ) =MU 

506  LET  AM  7 ,4 )  =  AM  o  ,3  j 

507  LET  AM(7,6)=3.0#MU 

bOo  LET  AM (.7 ,7)=AM(b,6) 

b09  LET  AM(7, 12)=2.0*MU 

510  LET  AMb,b)=AM7,4) 

b  11  LET  AM( 8 ,7 ) =2. 0#MU 

512  LET  AM(  8 ,  8)  =  AMI.  7,7) 

513  LET  ARK  8 » 1  b )  =  3 • 0#MU 

514  LET  AMI 9 , 8) =MU 

bib  LET  AM(y,9)=AM(8,b) 


B-lb 


bio 

517 

518 

519 

520 

521 

522 

523 

524 

525 
52b 
527 
52b 

529 

530 

531 

532 

533 

534 

535 

536 

537 
53b 

539 

540 

541 
54^ 

543 

544 

545 
54u 
547 
54b 

549 

550 

551 

552 

553 

554 

555 
557 
55b 
559 


FOR  N=4  TO  MLIM-1  DO 

LET  I=4#N-6  "AS  THE  ROW  INDEX 
LET  AMl,4*N-10)  =  (MP0P-N+1 )*LAMBDA 
LET  AM( I , I )=-( (MP0P-N)*LAMBDA+3.0*MU) 

LET  AM(I,4*N-1 )=MU 

ADD  1  TO  I  "I=4*N-5 

LET  AM  I  ,4*N-9) = (MPOP-N+1  )*LAMBDA 

LET  AM(I,I-1 )=3.0*MU 

LET  AM( I , I )  =  AM( 1-1 ,1-1 ) 

LET  AM(I,4*N)=2.0*MU 

ADD  1  TO  I  *'I=4*N-4 

LET  AM(l,4*N-b)= (MPOP-N+1 )*LAMBDA 

LET  AM( I , 1-1 ) =2 . 0*MU 

LET  AM( I , I ) =AM( 1-1 ,1-1 ) 

LET  AM( I , 4#N+ 1 ) =3 . 0*MU 

ADD  1  TO  I  "I  =  4*N-3 

LET  AM( I , 4*N-7 )  =  (MPOP-N+ 1 ) *LAMBDA 

LET  AM  1,1-1  )  =Mll 

LET  AM(I,I)=AM(I-1 ,1-1 ) 

LOOP  "OVER  THE  NUMBER  OF  CUSTOMERS  IN  THE  SISTEM 
' 'EQUATIONS  FOR  THE  LAST  FOUR  STATES. 

f  I 

LET  I-4«MLIM-6 

LET  AM(I,4*MLIM-10)=LAMBDA*(MPOP-MLIM+1 ) 

LET  AM( I , I ) =-3 .  0*MU 

ADD  1  TO  I  "I=4*MLIM-5 

LET  AM(l,4*MLIM-9)=LAMBDA*(MPOP-MLIM+1 ) 

LET  AM( 1,1-1 )=3.0*MU 
LET  AM(I,I)=AM(I-1 ,1-1) 

ADD  1  TO  I  "I=4*MLIM-4 

LET  AM(I,4*MLIM-S)=LAMBDA*UV1P0P-MLLM+1  ) 

LET  AMCI, 1-1 )=2.0*MU 
LET  AM( I , I) =AM( 1-1 , 1-1 ) 

ADD  1  TO  I  "Ir4*MLIM-3 

LET  AM(I ,4*MLIM-7 ) =LAMBDA*(MP0P-MLIM+1 ) 

LET  AM 1,1-1 ) =MU 

LET  AM( I , I ) =AM( 1-1 ,1-1 ) 

» 

1  OBTAIN  THE  INVERSE  OF  IN  PLACE. 

i 

'FACTOR  AM*,*). 


5b0 

561 

562 


564 

565 
5oo 
567 
56b 
5o9 
570 


CALL  SGEFA  (AM » ,* ) , IPVT(* ) , INFO) 

IF  INFO  NE  0 

PRINT  1  LINE  WITH  INFO  THUS 
TROUBLE  FACTORING  AM.  INFO  =  *. 

STOP 

OTHERWISE 

CALL  SG ED  I  (AM( » , *) , IPVT( * ) , DET( * ) , 1 1 ) 

| | OBTAIN  THE  SOLUTION  VECTOR  OF  THE  MATf IX  EQUATION. 
FOR  1=1  TO  MAXM,  LET  PV ( I ) =ML*AM( I , 1 ) 
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‘371 

372 

bib 

bin 

bib 

376 

377 
37a 

379 

380 
3a  1 
582 
383 
3a4 
333 
3au 
387 
533 
389 
340 

591 

bid 

bib 

394 

bio 

397 
598 
399 
000 
60  1 
602 
603 
b04 

605 

606 
□  07 
608 

609 

610 
Dll 
61 2 
613 
ol4 
b  15 
o  16 
6 1 1 
618 
o  19 
620 
621 

622 

623 


"SOLVE  FOR  THE  EMPTY  SYSTEM  STATE  PROBABILITY  (P.NULL)  AND  .CALCULATE 
' 'THE  CUSTOMER  STATE  PROBABILITY  VECTOR  (P.NSTATE(* ; ) . 

I  f 

LET  P.NSTATEC1  J  =  PV(1  )+PVl2)  "FOR  1  CUSTOMER 
LET  P . NSTATE (2)=PV(31+PV(4)+PV(5) 

LET  EN  SI S  =  P . NST  ATE ( 1 )+2.0»P.NSTATE(29 
LET  VNSYSrP. NSTATE ( 1 )+4 .0*P.NSTATE(2) 

LET  EN0=0.0 
LET  ESQ=0 . 0 

LET  SUM=P.NSTATE11 9+P.NSTATE129 
FOR  N  =  3  TO  MLIM  DO 
LET  I=4»N-6 

LET  P.NSTATE(N;=PV(I)+PV(I+1  )  +  PV(I+2)  +  PVU+3 ) 

ADD  P.NSTATE(N)  TO  SUM 
ADD  N*P.NSTATE(N;  TO  ENSYS 
ADD  (N-3)*P.NSTATE(N)  TO  ENw 
ADD  CN-39**2XP.NSTATE(N)  TO  ESw 
ADD  N**2*P.NSTATEtN)  TO  VNSYS 
LOOP  ' 'OVER  THE  OTHER  NON-ZERO  STATES 

f  I 

"CHECK  VALIDITY  OF  THE  PROBABILITY  SUM. 

I  t 

IF  SUM  LT  0.0  OR  SUM  GT  1.0 

PRINT  1  LINE  WITH  SUM  THUS 

ERROR  IN  ROUTINE  FINITE. ME2.Q.  PARTIAL  SUM  OF  STATE  PROBS  =  **.***». 
STOP 

OTHERWISE 

LET  P.NULL=1 .O-SUM 

LET  P.SYS.FULL  =  1  . O-P.NULL- P. NST ATE  ( 1  ) -P.  NSTATE  (2 )  "FOR  3  SERVERS 


t  t 
t  I 
I  » 


CALCULATE  THE  VARIANCE  AND  STANDARD  DEVIATION  OF  NUMBER  IN  THE  SYSTEM. 


LET  VNSYS= VNSYS- ENSYS* *2 
LET  SDNSYS=SQHT.F( VNSYS; 

LET  VNQ=ESQ-EN0**2 
LET  SDNG -SQRT . F ( VNQ ) 

I 

'CALCULATE  THE  MEAN  AND  STANDARD  DEVIATION  OF  BUSY  SERVERS. 

I 

LET  E.BUSY . SERVERS=P.NSTATE( 1 )+2.0*P.NSTATE(2 )+3.0*( 1 .O-P.NULL 
-P . NSTATE ( 1 ) -P . NSTATE (2 ) ; 

LET  VAR. BUSY. SERVERS=P.NSTATEl I )+4 .0*P.NSTATEC2 )+y .0* ( 1 .O-P.NULL 
-P .NSTATE ( 1 ) -P . NSTATE (2 ) ) -E . BUSY . SER  VERS**2 
LET  SD. BUSY. SER VERS=SQRT.FC VAR. BUSY. SERVERSJ 

\ 

’CALCULATE  MEAN  WAITING  TIME  USING  LITTLE’S  FORMULA. 

I 

LET  ESYS. WAIT-ENSYS/LAMBDA/C  REAL. FtMPOP) -ENSYS) 

LET  EQ. WAIT  =ENQ/LAMBDA/( REAL. F(MPOP) -ENSYS) 

IF  MP0P=3 
RETURN 

OTHERWISE  "CALCULATE  AND  PRINT  CONDITIONAL  PROB  DISTRIBUTIONS 
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• 'CALCULATE  THE  PROBABILITY  THAT  AN  ARRIVAL  FINDS  THE  SYSTEM  IN  STATE  N. 

t  t 

LET  NORM. CONST =M POP* P. NULL 
FOR  N= 1  TO  MLIM-1  DO 

ADD  (MPOP-N)*P.NSTATE(N)  TO  NORM. CONST 
LOOP  "TO  CALCULATE  THE  NORMALIZATION  CONSTANT 
LET  Q. NULL =MPOP*P. NULL/ NORM. CONST 
LET  QV ( 1 )=(MPOP-1 )*PV(1 ) /NORM. CONST 
LET  QV(2)=(MP0P-1 )*PV(2) /NORM. CONST 
LET  QV(3)=(MPOP-2)*PV (3) /NORM. CONST 
LET  Q V( 4 ) = (MPOP-2) *  PV (4 ) /NORM. CONST 
LET  QV (5 )=(MPOP-2)*PV(5) /NORM. CONST 
LET  Q.NSTATEC1 )=QV(1 )+QV(2) 

LET  Q.NSTATE (2 ) =QV C3 ) +QV (4 ) +QV (5) 

FOR  N=3  TO  MLIM-1  DO 
LET  I=4*N-6 

LET  y  V ( I )  =  ( MPOP- N ) *  PV ( I ) /NO RM . CONST 
LET  QV (1+1 )=(MPOP-N)*PV( 1+1 ) /NORM. CONST 
LET  Q V( 1+2 )  = (MPOP-N) *PV ( I+2J /NORM. CONST 
LET  QV C 1+3 )=CMP0P-N)*PV( 1+3) /NORM. CONST 
LET  Q. NSTATE(N)=QV(I)+QV(I+1 )+QV(I+2)+QV(I+3) 

LOOP  ' 'OVER  NUMBER  OF  ARRIVING  CUSTOMERS 

LET  P.  CUST .  WAIT=  1 .0-Q.NULL-Q.  NSTATEC 1  ) -Q.  NSTATEl.2 ) 

GO  TO  L2 

END  "ROUTINE  FINITE. ME2.Q 
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ROUTINE  SGEFA  (  A,  IPVT,  INFO; 


'’ROUTINE  FACTORS  THE  MATRIX  AC*,*)  INTO  UPPER  (U)  AND  STRICTLY  LOWER  (L) 
"TRIANGULAR  MATRICES  SUCH  THAT  A(*,»)  =  U(* , * )L( * , *  ) . 

''ROUTINE  IS  INTENDED  FOR  USE  WITH  OTHER  ROUTINES  OF  THE  LINEAR  OPERATIONS 
"  PACKAGE — LIN PACK.  THIS  VERSION  IS  A  CONVERSION  OF  THE  FORTRAN  ROUTINE 


''WRITTEN  BY  ClEVE  HOLER,  U.  OF  N.M.  AND 

»  I 

' ' ARGUMENTS: 

''NAME  MODE  ENTRY  VALUE 


''A  REALIN,  N)  SQUARE  MATRIX. 


' 'N  INTEGER  ORDER  OF  THE  MATRIX 

''IPVT  INTEGER(N) . 

''INFO  INTEGER  INDICATOR. 


ARGONNE  NAT  LAB. 


RETURN  VALUE 


AN  UPPER  TRIANGULAR  MATRIX  AND 
THE  MULTIPLIERS  WHICH  WERE  USED  TO 
GET  IT.  THESE  ARE  STORED  IN  L. 

A.  DIMENSION  OF  A(*,*). 

VECTOR  OF  PIVOT  INDICES. 

=  0  FOR  NORMAL  VALUE. 

=  K  IF  U(K,K;  EQ  0.0.  THIS 
INDICATES  THAT  SGESL  OR  SGEDI 
WILL  DIVIDE  BY  ZERO  IF  CALLED. 


DEFINE  I, INFO, J , K, KP1 , L ,N , NM1  AS  INTEGER  VARIABLES 
DEFINE  IPVT  AS  AN  INTEGER,  1 -DIMENSIONAL  ARRAY 
DEFINE  A  AS  A  REAL,  2-DIMENSIONAL  ARRAY 

I  » 

' 'GAUSSIAN  ELIMINATION  WITH  PARTIAL  PIVOTING. 

I  I 

LET  N=DIM.F(IPVT(*)) 

LET  INF0=0 
LET  NM1 =N-1 
IF  NM1  LT  1 
GO  TO  L7 
OTHERWISE 
FOR  .<  =  1  TO  NM1  DO 
LET  KP1  =K  +  1 


''FIND  L  =  PIVOT  INDEX  IN  THIS  COLUMN. 

I  » 

LET  SMAX=ABS.F(A(K,  K) ) 

LET  L=K 

FOR  I=K+1  TO  N  DO 

IF  ABS.F(A(I,Kj;  GT  SMAX 
LET  L=I 

LET  SMAXrABS.FU(I,X;  ; 

ALWAYS 

LOOP  ' ' FOR  MAX  ELEMENT 

LET  ipvt(k;=l 

I  I 

''ZERO  PIVOT  IMPLIES  THIS  COLUMN  ALREADY  TRIANGULAR I ZED. 

I  I 

if  a(l,k;  ^  o.o 

GO  TO  L4 
OTHERWISE 

i  i 

' ' INTERCHANGE  if  necessary. 
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IF  L  =  K 

GO  TO  LI 
OTHERWISE 
LET  T=A(L,K) 

LET  AIL,K1=A(K,K) 

LET  A(K,K)=T 

'Ll'  LET  T=-1.0/A(K,K) 

FOR  I =K  + 1  TO  N,  LET  A(I , K) =T*A( I ,K) 

(  I 

''ROW  ELIMINATION  WITH  COLUMN  INDEXING. 

I  ( 

FOR  J=KP1  TO  N  DO 
LET  T=A(L,J) 

IF  L=K 

GO  TO  L2 
OTHERWISE 
LET  A(L,J)=A(K,J) 

LET  AIK, J)=T 

'  L2'  FOR  I=K+1  TO  N,  LET  A(I,J)=T»AiI,K)+A(I,J) 

LOOP  "OVER  (J)  COLUMNS 
GO  TO  L5 

'L4'  LET  INFO=K 
' L5' LOOP  ' 'OVER  K 
’ L7 ' LET  IPVT(N)=N 
IF  A(N,N)=0.0 
LET  INFO=N 
ALWAYS 
RETURN 
END  ' ' SGEFA 
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ROUTINE  SGEDI  U,  IPVT,  DET ,  JOB) 


1  *  SGEDI  COMPUTES  THE  DETERMINANT  AND  INVERSE  OF  A  MATRIX  USING  THE  RESULTS 
' 'PRODUCED  BY  SGEFA. 

»  t 

’ ' ARGUMENTS: 

'  '  At  * , *  )  __  THE  REAL  FACTORED  MATRIX  FROM  SGEFA  ON  INPUT.  ON  OUTPUT  THE 
"  ARRAY  CONTAINS  THE  MATRIX  INV,  IF  REQUESTED.  OTHERWISE  UNCHANGED. 

' 1 IPVT  )  THE  INTEGER  PIVOT  VECTOR  FROM  SGEFA. 

' ’JOB  _  AN  INTEGER  SWITCH. 

"  ~  =  11  FOR  BOTH  DETERMINANT  AND  INVERSE. 

”  =01  FOR  INVERSE  ONLY. 

"  =10  FOR  DETERMINANT  ONLY. 

' 'DET(*>  _  CONTAINS  THE  DETERMINANT  OF  THE  MATRIX,  IF  REQUESTED.  OTHERWISE 
"  IS  NOT  REFERENCED.  THE  DETERMINANT  =  DET( 1 )* 10 .0**DET(2 ) ,  WITH 

' '  DET t 1 )  BETWEEN  0  AND  10,  AND  WITH  DET (2)  A  FLOATED  INTEGER. 

I  I 

’ 'NOTE:  A  DIVISION  BY  ZERO  WILL  OCCUR  IF  THE  INPUT  FACTOR  CONTAINS  A 

’’ZERO  ON  THE  DIAGONAL  AND  THE  INVERSE  IS  REQUESTED. 

I  I 

DEFINE  I,J,J0B,K,KB,KP1 ,L,N,NM1  AS  INTEGER  VARIABLES 
DEFINE  IPVT  AS  AN  INTEGER,  1 -DIMENSIONAL  ARRAY 
DEFINE  DET,  WORK  AS  REAL,  1 -DIMENSIONAL  ARRAYS 
DEFINE  A  AS  A  REAL,  2-DIMENSIONAL  ARRAY 
LET  N=DIM.F(IPVT(* ) ) 

RESERVE  WORKC*  >  AS  N  ’’ LOCALLY 

1  » 

' 'CALCULATE  THE  DETERMINANT  IF  REQUESTED. 

I  I 

IF  DIV.FtJOB,  10)  =  0 
GO  TO  Lu 
OTHERWISE 
LET  DET( I )=1 .0 
LET  DETC2 J  =0 .0 
LET  TEN=10.0 
FOR  1=1  TO  N  DO 

IF  IPVTU>  NE  I 

LET  DET( 1 )  =  -DET(1) 

ALWAYS 

LET  DETU  >=A(l,l)*DETU  > 

IF  DETU  >=0.0 
GO  TO  Lo 
OTHERWISE 

'Ll  '  IF  ABS.FCDETU  >>  GE  1.0 
GO  TO  L2 
OTHERWISE 

LET  DETU  >=TEN*DET(  1 ) 

SUBTRACT  1.0  FROM  DET(2> 

GO  TO  LI 

'  L2'  IF  ABS.FtDETU))  lt  ten 
GO  TO  L4 
OTHERWISE 

LET  DETU  )=DETU  >  /TEN 
ADD  1.0  TO  DETA2) 

GO  TO  L2 

'  LW '  LOOP  "OVER  I 
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'  'GET  INVERSE  OF  UPPER  TRIANGULAR  MATRIX  U(*,»). 

i  t 

' Lo ' IF  MOD. F( JOB,  101  =  0 
RELEASE  WORKC*) 

RETURN 
OTHERWISE 
FOR  K=1  TO  N  DO 

LET  A(K, K) =1  . 0/Al.K,  K) 

LET  T=-A(K,iO 

FOR  1=1  TO  K-1,  LET  A(I,K)=T*A(I,K) 

LET  KP1=K+1 
IF  N  LT  KPl 
GO  TO  L9 
OTHERWISE 
FOR  J=KP1  TO  N  DO 
LET  T=AlK,J) 

LET  A(K,J)=0.0 

FOR  1  =  1  TO  K,  LET  A(  I ,  J )  =T*A II , K)  +AQ ,  J ) 
LOOP  ' 'OVER  J 
' L9' LOOP  ' 'OVER  K 

t  t 

''FORM  INVERSE  ( Ul  *  INVERSE  (L.) 


LET  NM1 =N-1 
IF  NM1  LT  1 

RELEASE  WORK( * ) 

RETURN 

OTHERWISE 

FOR  KB= 1  TO  NM1  DO 
LET  K=N-KB 
LET  KPl =K+1 
FOR  I=KP1  TO  N  DO 

LET  WORK( I ) =A( I , K) 

LET  A( I , K) =0.0 
LOOP  ' 'OVER  I 
FOR  J=KP1  TO  N  DO 
LET  T=WORK( J) 

FOR  1=1  TO  N,  LET  A(I,K} =T*AU, J)+AII ,K) 
LOOP  ' 'OVER  J 
LET  L=IPVT(K1 

IF  l  NE  K  ''SWAP  ELEMENTS  OF  VECTORS  K  AND  L 
FOR  1=1  TO  N  DO 
LET  T=A(I,K) 

LET  A(I,K)=A(I,D 
LET  A(I,L)=T 
LOOP  "OVER  I  TO  SWAP 
ALWAYS 

LOOP  ' 'OVER  KB 
RELEASE  WORK(») 

RETURN 
END  "SGEDI 
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1  ROUTINE  MINRV  ( I PRINT, N , XRANGE, APARM, BPARM)  YIELDING  EMINX.VMINX 

2  '  ' 

j  " ROUTINE  CALCULATES  THE  MEAN  AND  VARIANCE  OF  THE  MINIMUM  OF  N  CONTINUOUS, 

4  ''POSITIVE,  I.I.D.  RANDOM  VARIABLES.  THE  C.D.F.  FOR  THESE  R.V.'S  IS 
B  ''CALLED  "FUN. CDF"  HERE,  AND  IS  SPECIFIED  BY  A  DEFINE- TO-MEAN  STATEMENT 
o  "IN  THE  PREAMBLE.  IF  THE  PRINT  FLAG  IPRINT  =1,  THE  BASIS  C.D.F.  AND 
7  "C.D.F.  OF  THE  MINIMUM  ARE  PRINTED  FROM  THIS  ROUTINE, 

b  '  ' 

9  "INPUTS: 

10  "IPRINT _ INTEGER  FLAG  TO  PRINT  FROM  THE  ROUTINE  (  =  1). 

11  "N _ NUMBER  OF  I.I.D.  RANDOM  IX-  )  VARIABLES  IN  THE  SET. 

IE  ' 'aRANGE  UPPER  LIMIT  ON  RANGE  OF  X,  USED  TO  CALCULATE  INTEGRATION  STEP. 
IS  "APARM  _ 21  FIRST  (REAL-VALUED)  PARAMETER  OF  THE  C.D.F.  OF  X. 

14  ' ' BP ARM  _  2  ND  ( REAL- VALUED;  PARAMETER  OF  THE  C.D.F.  OF  X. 

15  " 

10  "OUTPUTS: 

17  ' ' EMINX  _  EXPECTED  VALUE  OF  THE  MIN  OF  N  RANDOM  X-VARIABLES. 

1b  "VMINX _ 2_  VARIANCE  OF  THE  MIN  X. 

19  " 

20  DEFINE  I, IPRINT, J,K,M,N  AS  INTEGER  VARIABLES 

21  DEFINE  FV.SUMV  AS  REAL,  1 -DIMENSIONAL  ARRAYS 

22  RESERVE  FV ( * ) , SUMV ( * )  AS  2 

23  LET  il=2000  "STEPS  TO  INTEGRATE  ALLOWED 

24  LET  DELX =X RANGE/M  "INTEGRATION  STEP  SIZE  (SIMPSON'S  RULE) 

2b  IF  IPRINT  NE  1 

26  GO  TO  LO 

27  OTHERWISE 

26  SXIP  2  LINES 

29  PRINT  7  LINES  WITH  N, APARM, BPARM 

30  THUS 

PR03  DISTRIBUTION  FOR  THE  SMALLEST  OF  **  CONTINUOUS  RANDOM  VARIABLES 
PARAMETERS  OF  BASIS  DISTRIBUTION:  1ST  =  .  2ND  =  . 


ARGU¬ 

MENT 

BASIS  CDF  OF 

C.D.F.  MIN  SET 

3b 

' LO' LET 

FV(1  )  =  1  .0 

39 

LET 

FV(2)=0.0 

40 

FOR 

J  =  1  TO  2,  LET  3UMV( J )=FV(J) 

41 

FOR 

1=1  TO  M  DO 

42 

LET  X=I*DELX 

43 

IF  MOD. F ( I ,2 ) =0 

44 

LET  COEF  =2.0 

4b 

OTHERWISE 

4b 

LET  COEF  =4 . 0 

47 

ALWAYS 

4b 

LET  FX =FUW. CDF( APARM, BPARM,  X) 

4y 

LET  GN.C0MPL=(1 .0-FX)**N 

bO 

LET  FV(1 )=GN.COMPL 

b  1 

LET  FV(2 ) =X #GN . COMPL 

b2 

FOR  J=1  TO  2,  ADD  COEF*FV(J)  TO  SUMV(J) 

b3 

IF  HOD.FCI ,20)  =  U  AND  IPRINT=1  "PRINT  RESULTS 

o4 

PRINT  1  LINE  WITH  X , Fa , 1 . O-GN. COMPL 

bb 

THUS 

#.**»**  *,***** 
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37  ALWAYS 

58  IF  GN.COi'iPL  LT  0.0001  AND  COEF  =  2.0 

55  FOR  J=1  TO  2,  SUBTRACT  FV(J)  FROM  SUMV(J) 

60  GO  TO  LI 

o 1  OTHERWISE 

62  LOOP  "OVER  I 

63  ’Ll 'IF  IPRINT=1 

6M  PRINT  2  LINES  THUS 


67  ALWAYS 

68  LET  EMINX=SUMV(1 )*DELX/3.0 

65  LET  VMINX=2.0*SUMV(2)*DELX/5.0-EMINX**2 


70 

IF  I PRINT*  1 

71 

PRINT  2  LINES 

Id 

THUS 

MEAN 

VALUE  OF  MIN  R.V. 

STANDARD  DEV  MIN  R.V. 

75 

ALWAYS 

76 

RELEASE  FV(*) 

77 

RELEASE  SUMV(* ) 

78 

RETURN 

75 

END 

' 'MINRV 

WITH  EMINX,SQRT.F(VMIiiX) 

IN  THE  SET  _  . 

IN  THE  SET  . 
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ROUTINE  TO  LIMSTATE  ''OF  A  FINITE,  MULTI-SERVER  QUEUEING  SYSTEM’’  GIVEN 
LAMBDA,  MU,  MPOP,  NSERVE,  PLIM,  I  PRINT  YIELDING  LSTATE,  P.NULL,  P.NSTATE 

t  » 

' 'ROUTINE  SOLVES  FOR  THE  PROBABILITY  VECTOR  WHICH  CHARACTERIZES  THE  STEADY - 
' 'STATE  OF  A  QUEUEING  SYSTEM  WITH  A  FINITE  POPULATION  (MPOP)  OF  CUSTOMERS 
''SERVED  BY  NSERVE  CHANNELS  OF  EXPONENTIAL  SERVICE.  THE  ARRIVAL  RATE  FOR 
''EACH  INDIVIDUAL  NOT  IN  THE  SYSTEM  IS  A  CONSTANT  (LAMBDA).  THE  SERVICE 
'  'RATc,  IS  THE  CONSTANT  MU.  THE  ROUTINE  RETURNS  THE  (INTEGER)  STATE  INDEX 
' 'WHICH  IS  EXCEEDED  WITH  PROBABILITY  PLIM.  THE  ROUTINE  ALSO  RETURNS 
''THE  PROBABILITY  OF  AN  EMPTY  SERVICE  SYSTEM  (P.NULL)  AND  THE  VECTOR  OF 
''STATE  PROBABILITIES  (P .NSTATE(# ) ) .  REF:  GROSS  AND  HARRIS,  FUND.  QUEUE. 


' ' INPUT: 

’ ' LAMBDA 
'  'MU 
' 'MPOP 
' ' NSERVE 
' ’ PLIM 


' 'IPRINT 

I  I 

'  'OUTPUT: 

' ' LSTATE 

I  I 

' ' P.NULL 
' ' P.NSTATE 


I  I 


ARRIVAL  RATE  FOR  SERVICE  PER  INDIVIDUAL. 

SERVICE  RATE  PER  SERVER. 

CUSTOMER  POPULATION  SIZE. 

NUMBER  OF  SERVICE  CHANNELS. 

LIMIT  ON  THE  UPPER-TAIL  PROBABILITY  OF  THE  STATE 
PROBABILITY  DISTRIBUTION. 

INTEGER  SWITCH  TO  PRINT  FROM  THE  SUBROUTINE.  PRINTING 
OCCURS  WHEN  IPRINT=1. 

INDEX  OF  THE  MARKOV  STATE  WHICH  IS  EXCEEDED  WITH 
PROBABILITY  PLIM. 

PROBABILITY  THE  SERVICE  SYSTEM  IS  EMPTY. 

VECTOR  IN  WHICH  THE  N  TH  ELEMENT  IS  THE  PROBABILITY 
THAT  N  CUSTOMERS  ARE  IN  THE  SERVICE  SYSTEM. 


''ENDOGENOUS  VARIABLES: 

' ' ENO. DOWN  AVERAGE  NUMBER  OF  INDIVIDUALS  IN  THE  SERVICE  SYSTEM. 

' 'SD NO. DOWN  STD  DEV  NUMBER  OF  INDIVIDUALS  IN  THE  SERVICE  SYSTEM. 

’ ' P- SYS. FULL  PROBABILITY  THAT  ALL  SERVICE  CHANNELS  ARE  BUSY. 

"  E. BUSY. SERVERS  AVERAGE  NUMBER  OF  BUSY  SERVICE  CHANNELS. 

' ' SD. BUSY. SERVERS  STD  DEV  NUMBER  OF  BUSY  SERVICE  CHANNELS. 

' 'ENO. QUEUED  AVERAGE  NUMBER  OF  UNITS  WAITING  IN  THE  SERVICE  QUEUE. 

' 1 ESYS. WAIT  AVERAGE  WAITING  TIME  IN  THE  SERVICE  SYSTEM. 

'  '  EQ.  WAIT  AVERAGE.  WAITING  TIME  IN  THE  SERVICE  QUEUE,  INCLUDING 

' '  THE  INSTANCES  OF  ZERO  QUEUE  TIME. 


DEFINE  I, IPRINT, J, LSTATE, M, MAX, MPOP, N, NSERVE  AS  INTEGER  VARIABLES 
DEFINE  P.NSTATE  AS  A  REAL,  1 -DIMENSIONAL  ARRAY 
RESERVE  P.NSTATE(»)  AS  MPOP 
LET  R=LAMBDA/MU 
LET  CONSERVE 
IF  NSERVE  GE  MPOP 

PRINT  1  LINE  W ITH  NSERVE, MPOP 
THUS 

INPUT  ERROR  TO  ROUTINE  LIMSTATE.  NSERVE  =  »»#*»  MPOP  =  »»»»». 
RETURN 
OTHERWISE 

LET  CDF. LIM= 1.0- PLIM 

LET  RHO=R/C 

LET  ENO. DOWN =0.0 

LET  E.BUSY .SERVERS=0.0 

LET  VAR. BUSY. SERVERS=0.0 

LET  SUM= 1 . 0 
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LET  QSUM=C 
LET  NFACTORIAL=1 . 0 
LET  RATIO. FACTORIALS^  .0 
IF  NSERVE=1 

GO  TO  PASS 
OTHERWISE 

FOR  N=1  TO  NSERVE-1  DO 

LET  NFACTORIAL=NFACTORIAL*N 

LET  RATIO. FACT0RIALS=RATI0.FACT0RIALS«(MP0P-N+1 ) 

LET  P.N=R"*N/NFACTORIAL*RATIO. FACTORIALS  "OMITTING  "P.NULL 
LET  P.  NSTATECN )  =  P. N  "WITH  THE  NTH  STATE  STORED  IN  THE  NTH  ELEMENT 
ADD  N*P.  N  TO  ENO . DOWN 
ADD  (nserve-n;*p.n  TO  WSUM 
ADD  P.N  TO  SUM 
ADD  N*P.N  TO  E. BUSY. SERVERS 
ADD  N**2"P.N  TO  VAR . BUSY . SERVERS 
LOOP  ' ' TO  DEVELOP  SUMS 

t  I 

"CALCULATE  THE  LAST  TERM  IN  THE  SUM  FOR  P.NULL 

t  I 

' PASS' LET  CFACTORI AL=C"NFACTORIAL 
LET  CEXPOC=C*"N SERVE 
LET  COEF  =CEXPOC/C FACTORIAL 

LET  SUM. LAST  =0.0  "SUM  FOR  PROB  THAT  NO  IN  SYSTEM  GT  OR  =  NSERVE 
LET  VNO.DOWN=VAR. BUSY. SERVERS 
FOR  N=NSERVE  TO  MPOP  DO 

LET  RATIO . F ACTORI ALS  =  RATIO . F ACTORI ALS" { MPOP-N+ 1 ) 

LET  P. N=COEF*RATIO.FACTORIALS"RHO"*N 
ADD  P.N  TO  SUM. LAST  "OMITTING  "P.NULL 
ADD  N*P.  N  TO  ENO. DOWN  "OMITTING  "P.NULL 
ADD  N«*2«P.N  TO  VNO.DOWN  "OMITTING  "P.NULL 
LET  P.NSTATE(N)=P.N  "OMITTING  "P.NULL 
LOOP  "TO  OBTAIN  THE  LAST  TERM 
ADD  SUM. LAST  TO  SUM 

LET  P.NULL=1 . 0/SUM  "PROB  OF  SYSTEM  NULL  STATE 

f  t 

’’CHECK  OF  CALCULATED  SERVICE  SYSTEM  STATE-PROBABILITIES. 

i  i 

IF  P.NULL  GE  1.0  OR  P.NULL  LE  0.0 
PRINT  1  LINE  WITH  P.NULL  THUS 

ERROR  IN  CALCULATING  STATE  PROBABILITIES.  P.NULL  =  . 

STOP 

OTHERWISE  "CALCULATE  EXPECTED  VALUES 
LET  ENO. DOWN=EHO.DOWN«P. NULL 
LET  VNO . DOWN=VNO . DOWN«P . NULL- ENO . D0WN«*2 
LET  SD NO . DOWN  =S0 RT . F ( VNO . DOWN ; 

LET  P. SYS. FULL=P.NULL«SUM. LAST 

LET  E. BUSY. SERVERS=P.NULL«E. BUSY. SERVERS+C«P. SYS. FULL 

LET  VAR. BUSY. SERVERS=P. NULL*VAR . BUSY. SERVERS+C«*2*P. SYS. FULL 

-E . BUSY . SERVERS*"^ 

LET  SD . BUSY . SERVERS=SURT. F (. VAR . BUSY . SERVERS; 

LET  ENO. QUEUED=ENO. DOWN-NSERVE+P. NULL *0 SUM 

LET  ESYS.WAITi:ENO.  DOWN/LAMBDA/ (  REAL. F(.MPOP) -ENO.  DOWN; 

"ABOVE  EXPRESSION  IS  KNOWN  AS  LITTLE'S  FORMULA. 

LET  EO.WAIT=ESYS. WAIT-1 .O/MU 
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112 

113 

114 
11 ‘3 

1  1o 
117 
1  18 
11* 

120 
121 
12<i 
1 23 
12A 
12;) 

12o 
127 

STATISTICS  FOR  THE  STEADY  STATE  OF  THE  SERVICE  SYSTEM 
130  PRINT  o  EINES  WITH  ENO. DOWN, SDNO. DOWN, P. SYS. FULL, E. BUSY. SERVERS, 


131  3D. BUSY. SERVERS, ENO. QUEU ED, ESYS. WAIT, £Q. WAIT 

132  THUS 

AVERAGE  NUMBER  OF  INDIVIDUALS  IN  THE  SERVICE  SYSTEM  _  **#.*#** 

STD  DEV  NUMBER  OF  INDIVIDUALS  IN  THE  SERVICE  SYSTEM  _  ##*.###* 

PROBABILITY  THAT  ALL  SERVICE  CHANNELS  ARE  BUSY  *.*#*# 

AVERAGE  NUMBER  OF  BUSY  SERVICE  CHANNELS  *.*#*# 

STD  DEV  NUMBER  OF  BUS7  SERVICE  CHANNELS  ».*#** 

AVERAGE  NUMBER  OF  INDIVIDUALS  QUEUED  FOR  SERVICE  *##.*### 

AVERAGE  WAITING  TIME  IN  THE  SERVICE  SYSTEM  #»##.### 

AVERAGE  WAITING  TIME  IN  THE  SERVICE  QUEUE  ****.*#» 

141  SKIP  2  LINES 

14^  RETURN 

143  END  "LIMSTATE 


''CALCULATE  SERVICE  SYSTEM  STATE  PROBABILITIES. 

I  I 

LET  P3UM=P . NULL 
FOR  N= 1  TO  MPOP  DO 

LET  P.NSTATEt,N;  =  P.NSTATEU\l)*P.NULL 
ADD  P.NSTATE(N)  TO  PSUM 
IF  PSUM  LE  CDF.LIM 
LET  LSTAT£=N 
ALWAYS 

LOOP  ’’OVER  ALL  STATES 
IF  IPRINT  WE  1 
RETURN 
OTHERWISE 
SKIP  2  LINES 
PRINT  2  LINES  THUS 


t 
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